REPORT  DOCUMENTAT 


AFRL-SR-BL-TR-98-- 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  a 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewint 
collection  of  information,  including  suggestions  for  reducing  this  burden, 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office 


1.  AGENCY  USE  ONLY  (Leave  blank) 

2.  REPORT  DATE 

18  June  1998 

3.  REPORT  TYPE  AND  DATES  COVERED 

Final  Technical  Report  01  May  94  to  30  Apr  98 

4.  TITLE  AND  SUBTITLE 

Fracture  and  Failure  at  and  Near  Interfaces  Under  Pressure 

5.  FUNDING  NUMBERS 

F49620-94- 1-0253 

6.  AUTHOR(S) 

W.  G.  Knauss,  PI 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Graduate  Aeronautical  Laboratories 

California  Institute  of  Technolocy 

Pasadena,  CA  91125 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR/NA 

110  Duncan  Avenue,  Ste  B115 

Bolling  AFB,  DC  20332-8050 

10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

F49620-94-1-0253 

11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 
Approved  for  public  release;  distribution  unlimited 


13.  ABSTRACT  (Maximum  200  words) 

This  work  addressed  the  failure  behavior  of  solid  propellant  rocket  fuels  through  crack  propagation.  The  objective  of  the 
study  was  to  1)  develop  the  means  for  measuring  large  deformation  fields  around  the  tips  of  stationary  or  slowly  moving 
cracks,  to  develop  realistic  data  for  comparison  with  improved  analytical  results,  and  to  2)  initiate  a  new  computational 
approach  for  stress  analysis  of  cracks  at  and  near  interfaces,  which  can  draw  on  the  expanding  capabilities  of  parallel 
processing. 

Important  results  are  that 

1)  strain  inhomogeneities  are  much  more  pronounced  than  hitherto  anticipated;  they  are  associated  with  the  granular 
microstructure  and  are  characterized  by  spatial  variations  on  the  sub-millimeter  size  scale; 

2)  these  strain  inhomogenieties  dominate  the  deformation  field  around  a  crack  tip  and  control  the  fracture  process; 
cracks  connect  statistically  distributed  regions  of  high  strain,  if  they  are  reasonably  close  to  the  crack  propagation  line; 

3)  a  domain  of  1-3  mm  in  extent  around  the  crack  tip  is  totally  dominated  by  this  highly  inhomogeniously  deformed 
material,  with  the  adjacent  material  undergoing  large  nonlinear  deformations; 

4)  the  disintegrating  region  close  to  the  crack  tip  is  not  likely  to  be  describable  in  continuum  terms,  and  needs  to  be 
modeled  in  a  discrete  fashion. 


14.  SUBJECT  TERMS 

15.  NUMBER  OF  PAGES 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

18.  SECURITY  CLASSIFICATION 

19.  SECURITY  CLASSIFICATION 

20.  LIMITATION  OF 

OF  REPORT 

OF  THIS  PAGE 

OF  ABSTRACT 

ABSTRACT 

Unclassified 

Unclassified 

Unclassified 

UL 

ictions,  searching  existing  data  sources, 
den  estimate  or  any  other  aspect  of  this 
Operations  and  Reports,  1215  Jefferson 
-0188),  Washington,  DC  20503. 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


Final  Technical  Report 
to  the 

Airforce  Office  of  Scientific  Research 
on  Grant  No.  F49620-94-1-0253 


Fracture  and  Failure  at  and  Near 
Interfaces  Under  Pressure 


W.  G.  Knauss,  PI 

June  1998 

Graduate  Aeronautical  Laboratories 
California  Institute  of  Technology 

Pasadena,  California  91125 


PTIC  QUALITY  DJSPECTED  1 


Fracture  and  Failure  at  and  Near  Interfaces  Under  Pressure 


Abstract 

A  research  study  is  summarized  that  addressed  the  failure  behavior  of  solid  propellant 

rocket  fuels  through  crack  propagation.  The  objective  of  the  study  was  to 

1)  develop  the  means  for  measuring  large  deformation  fields  around  the  tips  of 
stationaiy  or  slowly  moving  cracks,  in  order  to  develop  realistic  data  against 
which  improved  analytical  results  can  be  compared,  and  to 

2)  initiate  a  new  computational  approach  for  stress  analysis  of  cracks  in  solid 
propellants  at  and  near  interfaces,  which  analysis  can  draw  on  the  ever 
expanding  capabilities  of  computationally  parallel  processing. 

Correspondingly,  both  experimental  and  analytical/numerical  work  was  carried  out: 

(A)  Experimental  work: 

(1)  A  full  field  method  for  visualizing  and  measuring  deforaiation  around  the  crack 
tip  in  a  fracture  process  with  large  stoins  has  been  developed; 

(2)  The  method  was  applied  to  the  fracture  in  a  solid  propellant  (growth  of  a  crack). 
It  was  found  that  lai'ge  variations  of  strain  occur  in  nominally  “homogeneous” 
deformation  fields,  with  sttain  variations  by  as  much  as  factors  of  three. 

(B)  Computational/Analytical  work: 

(1)  Proposed  and  implemented  a  discrete-continuum  approach  to  analyze  damage 
evolution  in  solid  propellant  (particle-filled  elastomers  in  general). 

(2)  Established  analytical  stability  conditions  for  interfacial  dewetting  at  planar, 
cylindrical  and  spherical  interfaces  for  uniforai  decohesion. 

(3)  Applied  the  method  to  investigate  the  effect  of  particle  interaction  on  the  damage 
evolution  in  solid  propellant:  Small  changes  in  particle  geometry  and  spacing 
have  large  effects  on  pai'ticle  debonding.  Reattachments  observed  as  part  of  the 
global  dewetting  process. 

An  important  result  of  this  study  is  that 

1)  strain  inhomogeneities  are  much  more  pronounced  than  hitherto  anticipated; 
they  are  associated  with  the  granular  microstnicture  of  the  material  and  are 
chai-acterized  by  spatial  variations  on  the  sub-millimeter  size  scale; 

2)  these  strain  inhomogeneities  dominate  the  deformation  field  around  a  crack  tip 
and  control  the  fracture  process,  in  that  cracks  tend  to  connect  statistically 
distributed  regions  of  high  strain,  if  they  ai'e  reasonably  close  to  the  crack 
propagation  line; 

3)  a  domain  of  1-3  mm  in  extent  around  the  crack  tip  is  totally  dominated  by  this 
highly  inhomogeniously  deformed  material,  with  the  adjacent  material 
undergoing  large  nonlinear  deformations; 

4)  the  disintegrating  region  close  to  the  crack  tip  is  not  likely  to  be  describable  in 
continuum  teims,  and  needs  to  be  modeled  in  a  discrete  fashion. 


1.  Introduction 


Since  the  development  of  solid  propellant  rocket  motors  in  the  late  1950’,  an  important 
operational  and  design  consideration  has  always  been  the  reliable  functioning  and  readiness 
of  both  tactical  and  sti-ategic  missile  systems.  The  most  impoitant  and  most  difficult 
component  of  the  system  analysis  has  been  the  predictability  or  suppression  of  failure  by 
the  intitation  and/or  propagation  of  cracks  in  the  fuel  charge  or  grain.  Such  failures  typically 
arise  in  the  grain  at  “star  points”  or  at  the  support  interface  of  the  grain  to  the  motor  case. 
Either  situation  is  equally  detiimental,  because  crack  initiation  can  lead  to  propagation  as  a 
result  of  the  highly  transient  pressure  pulse  during  ignition,  or  during  long  term  storage 
while  awaiting  deployment.  For  tactical  missiles  the  exposure  to  under-wing  use  during 
readiness  maneuvers  provides  an  additional  deliterious  use  environment. 

Two  aspects  of  fracture  analysis  and  prediction  in  composite  propellants  have  concerned 
the  mechanics  community  ever  since  the  start  of  solid  propellant  use,  and  for  neither  has  a 
satisfactory  solution  been  found  in  the  roughly  four  inteiwening  decades:  One  addresses  the 
experimental  deteiTnination  of  the  strain  fields  around  the  tips  of  cracks  in  these  materials; 
the  other  is  the  computation  of  stress  fields  either  form  these  strains  or  from  other  “first” 
principles.  Virtually  all  analyses  have  been  based  on  variants  of  linearly  elastic  analyses, 
which  have  long  been  known  to  provide  estimates  at  best,  sinee  the  nonlinear 
characteristics  of  the  propellant  materials  have  been  well  established  for  several  decades. 
However,  the  means  of  dealing  with  that  type  of  material  charaeterization  is  not  at  all  well 
in  hand  nor  has  it  been  developed  for  applications  to  motor  design. 

1.1  A  new  way  for  stress  analysis  involving  high  strain  (stress)  gradients. 

For  over  three  decades  propellant  stress  analysis  has  been  accomplished  by  means  of  finite 
element  codes  that  started  with  lineai'ly  elastic  material  chai'acterization.  While  specialty 
codes  have  been  developed,  to  address  nonlinear  behavior  of  the  material,  they  h'ave  neither 
been  demonsti-ated  in  connection  with  appropriate  experimental  input,  nor  are  the 
constitutive  equations  used  for  them  particularly  gennane  to  propellants.  All  these 
computations  have  to  be  supplemented  with  a  heavy  dose  of  empiricism  based  on  years  of 
practical  experience  to  make  failure  estimation  roughly  comparable  with  motor 
performance. 


The  most  common  assumption  underlying  these  analyses  is  that,  in  spite  of  its 
microstructure,  the  propellant  is  an  essentially  homogeneous  material.  The  consequence  of 
this  assumption  is  that  strain  fields  are  smoothly  varying  and  maxima  —  drivers  of  failure 
initiation  —  occur  where  continuum  theory  says  they  should  occur.  We  shall  see  that  the 
experimental  analysis  of  this  question  tells  us  otherwise,  and  the  variations  in  strain  fields 
arising  from  inhomogeneities  is  indeed  very  significant.  The  consequence  for  failure 
prediction  is  then  that  such  behavior  needs  to  be  incoiporated  into  analyses.  The  question 
of  how  to  do  this  is  not  answered  simply. 

Typically  propellants  have  been  modeled  as  (homogeneous)  continua.  Tliat  approximation 
seem  reasonable  when  the  strain  variations  are  not  important  in  the  fracture  sense  and  when 
the  domain  considered  is  large  in  comparison  to  the  domains  in  which  inhomogeneities 
dominate.  However,  when  these  latter  domain  sizes  are  on  the  order  of  domains  that 
determine  fracture,  such  as  at  the  front  of  a  crack,  then  continuum  representations  are  no 
longer  appropriate.  Attempts  have  been  made  over  the  past  years  to  characterize  the 
disintegrating  material  at  tlie  crack  tip  as  a  continuum  with  highly  nonlinear  properties,  thus 
the  odd  situation  arises  whereby  one  discretizes  the  crack  tip  domain,  tries  to  find  a 
“smeared-out”  ore  average  continuum  model  of  the  failing  material,  and  then  derives  from 
that  a  failure  criterion  that  comes  close  to  the  physical  situation  that  is  timly  discreet  in  its 
failure  aspects. 

Nonlinear  behavior  is  widely  observed  in  particulate  composites,  such  as  solid  propellants, 
tires,  toughened  plastics,  even  when  global  defoitnations  are  relatively  small.  Global 
nonlinearity  in  particulate  composites  is  often  caused  by  cavitation  or  tear  (appearance  of 
voids  in  the  matrix)  and/or  by  interfacial  particle  debonding  or  "dewetting";  see 
Schippel(1920),  Smith(1959),  Farris(1964),  Oberth  (1967),  Knauss  and  Mueller  (1979) 
and  Gent  and  Pai-k  (1984).  Dewetting  and  interfacial  void  generation  occurs  when  the  bond 
at  an  interface  is  relatively  weak;  otherwise  cavitation  can  occur.  How  to  represent  this 
damaged-material  behavior  in  continuum  terms  for  finite  element  analyses  has  been  a  long 
standing  question.  Constitutive  models  containing  damage  have  been  proposed,  for 
example,  by  Fairis  and  Schapeiy  (1973),  Schapery  (1986),  Govindjee  and  Simo  (1992), 
Vratsanos  and  Farris  (1993)  and  Ravichandran  and  Liu  (1995).  However,  these  models, 
excepting  that  Govindjee  and  Simo  (1992),  are  typically  restricted  to  infinitesimal 
elastic/viscoelastic  matrix  behavior,  and  therefore  their  applications  in  regions  of  high 
deformation  gradients  such  as  a  crack  tip,  are,  at  best,  questionable.  More  importantly, 
one  must  critically  examine  the  implied  proposition  of  first  representing  a  tmly 


discontinuous  process  by  a  homogenized  continuum  formulation  and  then  turn  around  and 
formulate  a  new  local  failure  criterion  that  addresses  the  crack  propagation  process  by 
means  of  the  homogenized  material  description  in  such  a  way  that  the  true  physical  process 
of  the  discrete  micro  crack  growth  and  coalescence  is  reproduced  faithfully  near  the  crack 
tip. 


Figure  1;  Above:  Tripartition  of  the  crack  tip  domain,  with  local  discrete  behavior  at  the 
crack  tip,  nonlinear  continuum  response  farther  away,  and  linear  behavior  over  the  rest  of 
the  large  solid.  Below:  Detailed  account  of  the  crack  tip  region  close  to  the  crack  tip.  Daric 
areas  represent  voids,  gray  areas  binder  and  white  particles. 


An  alternate  way  to  crack  tip  characterization,  that  circumvents  this  circuitous  route  for 
computing  involving  two  approximation  schemes  is,  thus,  to  feat  the  zone  around  the 
crack  tip  in  a  discrete  fashion,  that  is,  ultimately  with  statistical  particle  distributions,  and 
deal  with  the  failure  around  individual  pai'ticles. 

It  appears,  therefore,  more  reasonable  to  characterize  the  material  locally  near  a  crack  tip  in 
a  discrete  manner  by  modeling  a  select  region  of  the  high  strain  gradient  domain  through 
discretely  failing  elements  embedded  in  another  and  larger  region  of  nonlinear  and/or  linear 
material  response.  Here  failure  sites  ultimately  become  the  fracture  sites  by  which  a  crack 
propagates  through  void  coalescence.  Such  a  hybrid  discrete-continuum  approach  appears 
feasible  in  light  of  the  ever  more  rapidly  increasing  power  of  computing  machines.  We  thus 
pursue  here  initially  the  time-independent  (non-viscoelastic)  problem  of  the  dewetting 
process  allowing  for  lai-ge  defomiation  in  the  matrix  material,  and  with  tlie  intention  of 
incorporating  the  results  into  a  larger  finite  element  analysis  for  crack  growth  studies  in 
particulate  composites. 

Figure  1  shows  such  an  approach  graphically:  A  small  zone  around  the  crack  tip  is  modeled 
with  individual  particles,  the  next  farther  (but  “close  in”)  domain  from  the  crack  tip  as  a 
continuum  possessing  nonlineai-  characteristics,  and  the  remainder  of  the  stmcture  is 
sufficiently  well  characterized  by  linear  behavior  since  the  strains  tend  to  be  rather  small 
there. 


1,2.  Experimental  definition  of  the  discrete  analysis  domain 

Experimental  observations  are  necessary  to  assess  the  size  of  the  domain  within  which 
dicretization  of  particulate  mechanics  is  appropriate.  This  endeavor  was  the  first  pait  of  the 
project  and  revolved  around  microscope  studies  of  deformations  and  strain  determination 
through  the  Digital  Image  Correlation  Method  (DIG).  In  this  we  concentrated  on  the 
microdamage  evolution  and  its  interaction  with  small  cracks.  One  observes  clearly  that  the 
deformation  in  solid  propellants  is  inhomogeneous  at  the  length  scale  of  particles 
(especially  in  high  defonnation  regions,  e.g.  near  a  crack  tip),  while  tlie  defoitnation  is 
“homogeneous”  at  a  lai'ger  length  scale.  Experimental  observations  indicate  that 


microdamage-induced  heterogeneity  plays  the  dominant  role  in  crack  propagation  and  the 
strain  distribution  field  near  a  crack  tip. 


The  Digital  Image  Correlation  Method. 

Because  the  modifications  to  the  original  Digital  Image  Correlation  code  (DIC)  sought  here 
are  closely  tied  to  the  inner  workings  of  the  algorithm  for  the  two-dimensional  correlation 
scheme,  the  latter  is  summarized  briefly.  A  three-dimensional  extension  to  the  code  is 
presented  in  “Submicron  Defomiation  Field  Measurements  Part  II:  Improved  Digital  Image 
Correlation”  by  Vendroux,  G.  and  Knauss,  W.G  ,Exper mental  Mechanics,  38,  No.  2,  pp. 
86-92,  June  1998. 


Two-dimensional  Digital  Image  Correlation 

A  surface  profile,  as  obtained,  for  example,  by  a  Scanning  Tunneling  Microscope,  is  a 
discrete  record  of  the  "height"  of  the  surface  at  grid  points  assigned  to  a  specimen  surface. 
Let  f(x,y)  represent  the  surface  profile  of  a  specimen  in  an  undefonned  state  at  point 
G(x,y),  and  g(x,  y)  the  surface  profile  after  defomiation  at  the  corresponding  point 
G{x,  y).  If  the  profile  pattern  before  defonuation  is  uniquely  related  to  the  profile  pattern 
after  deformation,  a  coirelation  of  these  two  patterns  exists  to  detect  the  profile  difference 
which  represents  the  object  defonuation.  Let  %  be  the  mapping  from  the  undefonned  to  the 
deformed  configuration 

G->G=x(G)  such  that  g(x,  y)  =  f(x,y)  (1) 


x=x+u(x,y)  (2) 

y  =  y  +  v(x,y) 


with  u  and  v  the  in-plane  displacements  of  G,  and  let  Go(xo,5'o)  be  the  image  of  GoCxo.y,,) 
through  %;  further,  let  S  be  a  (sub)set  of  points  around  G^  and  S  the  corresponding  (sub)set 
of  point  around  Gg.  Assuming  that  S  is  sufficiently  small,  eq  (2)  can  be  expanded  into 


•X  —  X  +  H(Xo,yo)  ^  I  Xq  )  ^  I  X(, (y  3^0  ) 

X  =  y  +  v(.Xo,  yo) + £|  (x  -  Xo )  +  -^1  (y  -  yo ) 


du[ 


(3a) 

(3b) 


as  the  linearization  of  %  around  Gg.  For  a  discrete  set  of  data  define  the  correlation 
coefficients 


(4a) 


C  = 


2„,.,[/(G,)-s(i,(G,))r 


or 


C  =  l-- 


(4b) 


It  is  clear  that  C  will  be  zero  when  the  coefficients  of  the  mapping  %,{....}  are  indeed  the 
displacements  and  their  derivatives  at  Gq  [4,  8].  The  best  estimate  of  these  values  are 
determined  by  minimizing  C,  which  process  can  be  viewed  as  a  non-linear  optimization 
scheme,  some  details  of  which  will  be  discussed  in  section  3.2  under  “Optimization 
Scheme”. 


2.  Analytical 

The  numerical  work  has  centered  on  developing  segment  or  elements  containing  particles, 
which  undergo  the  dewetting  process.  The  issues  for  resolution  were  primarily  the  stability 
of  the  dewetting  process  in  the  computational  model.  We  believe  that  this  problem  is  now 
well  in  hand.  In  a  similar  vain,  the  stability  of  cracks  branching  from  the  dewetting  void  to 
fracture  the  binder  is  similarly  handled.  As  a  result  of  this  research  it  has  also  become  clear 
that  a  somewhat  coarser  modeling  of  the  dewetting/fractum  process  is  tolerable.  It  is,  by 
now,  abundantly  clear  that  tlie  proposed  scheme  of  highly  detailed  particle  induced  failure 
around  crack  tips  can  be  treated  numerically  with  the  currently  available  (parallel) 
computing  systems.  Moreover,  problems  smaller  than  full-scale  motor  applications  can  be 
addresses  probably  without  parallel  processing.  This  assessment  is  based  on  our 
experience  that  less  refined  damage  models  around  pai-ticles  can  now  be  constracted. 


3.  Report  content 


In  the  following  development  the  experimental  work  is  presented  first.  That  presentation  is 
followed  by  a  delineation  of  the  numerical  modeling.  However,  because  these  studies  have 
been  already  documented  in  written  reports  and  publications  they  are  only  abstracted  and 
summarized  here,  with  the  full  papers  denoted  by  an  asterisk  attached  as  appendices. 

3.1.  Strain  Inhomogeneitv  and  Discontinuous  Crack  Growth  in  a 
Particulate  Composite. 

This  topic  has  been  documented  as  an  Aeronautical  Engineer’s  Thesis  under  the  title: 
Gonzales,  J.; 

Full  Field  study  of  Strain  Distribution  near  the  Crack  Tip  in  the  Fracture  of 
Solid  Propellants  via  Large  Strain  Digital  Image  correlation  and  Optical 
Microscopy;  Engineer’s  Thesis,  1996,  California  Institute  of  Technology,  Pasadena,  CA 
91125 

Abstract 

A  full  field  method  for  visualizing  defomation  around  the  crack  tip  in  a  fracture  process 
with  large  strains  is  developed.  A  digital  image  coirelation  program  (DIC)  is  used  to 
incrementally  compute  stains  and  displacements  between  two  consecutive  images  of  a 
deformation  process.  Values  of  strain  and  displacements  for  consecutive  deformations  are 
added,  this  way  solving  convergence  problems  in  tlie  DIC  algorithm  when  lai-ge 
deformations  aie  investigated.  The  method  developed  is  used  to  investigate  the  strain 
distribution  within  1  mm  of  the  crack  tip  in  a  pai'ticulate  composite  solid  (propellant)  using 
microscopic  visualization  of  the  defoimation  process. 


A  paper  based  on  that  work  will  appear  in  September  1998  in  the  International  Journal  of 
Solids  and  Structures: 

Gonzales,  J.  and  Knauss,  W.G.; 

Strain  Inhomogeneity  and  Discontinuous  Crack  Growth  in  a  Particulate 
Composite*; 

Abstract 

A  full  field  method  for  visualizing  defoimation  around  tlie  crack  tip  in  a  fracture  process 
with  large  strains  is  developed.  Digital  Image  Coirelation  (DIC)  is  used  to  incrementally 
compute  strain  and  displacements  between  consecutive  images  in  a  global  deformation 
process.  Values  of  strains  and  displacements  for  consecutive  defoimations  are  added,  this 


way  solving  convergence  problems  in  DIC  algorithm  when  large  deformation  are 
investigated.  The  method  developed  here  is  used  to  investigate  the  strain  distribution  within 
1mm  of  the  crack  tip  in  a  particulate  composite  solid  propellant.  It  is  shown  that  in 
“nominally”  homogeneous  defonnations  (unidirectional  deformation  of  a  propellant  sheet) 
the  strains  vary  over  a  range  of  a  factor  of  three,  with  the  average  corresponding  to  the 
global  strain  defined  as  the  ration  of  boundary  displacement  and  the  specimen  width  or 
height.  In  crack  propagation  studies  such  high  strain  regions  are  the  loci  for  crack 
propagation. 

A  second  presentation  of  essentially  the  same  material  entitled 

Strain  Distribution  around  Cracks  in  Damaged  Particle-Filled  Elastomers 
has  been  made  by  J.  Gonzales  at  the  1997  Composite  Conference  held  at  Hawaii  in  July 
1997. 

The  standard  use  of  the  Digital  Image  Coixelation  (DIC)  works  well  to  investigate 
deformations  where  tlie  maximum  principal  strain  is  less  than  10%.  For  studying 
deformation  when  this  limit  is  reached,  it  is  necessary  to  use  additional  tools  like  the  Large 
Deformation  Digital  Image  Con-elation  (LD-DIC).  This  method  perfoi-ms  well  using  the 
standard  DIC  in  a  stepwise  manner,  following  the  defoi-mation  history.  Results  for  two 
tests  on  a  particle-filled  elastomer  are  presented  using  the  LD-DIC.  The  strain  associated 
with  void-opening  in  a  particle-filled  elastomer  (solid  propellant)  is  addressed  in  the  first 
test  while  the  effect  of  damage  on  the  strain  distribution  in  a  crack  opening  problem  for  the 
same  material  is  investigated  in  the  second  situation. 


3.2.  The  next  papers  cover  the  numerical/analvtical  work  on  the  deformation  of 
damage  mechanics  in  particulate  filled  composites. 

This  phase  of  the  research  addresses  tlie  micromechanics  of  a  dewetting  particulate 
composite  to  model  both  the  failure  progression  as  well  as  the  global  force-displacement 
(stress-strain)  response.  The  idea  is  that  such  a  model  becomes  the  small-scale  finite 
element”  for  a  larger  propellant  domain  in  which  damage  and  failure  progress  as  coupled  to 
the  nonlinear  response  of  the  local  material. 


The  first  paper  is 

Zhong,  X.A.  and  Knauss,  W.G. 

Analysis  of  Interfacial  Failure  in  Particle-Filled  Elastomers* 

Journal  of  Engineering  Materials  and  Technology  (ASME  Transactions);  119,  pp.  198- 
204,  July  1997  (presented  at  the  1997  ASME/ASCE/SES  Joint  Summer  Meeting  at 
Northwestern  University,  McNu'97) 

Abstract 

The  evolution  of  microdamage  (interfacial  dewetting)  in  highly  filled  elastomers  under 
consideration  of  high  defonnation  gradients  is  examined  with  the  aid  of  the  ABAQUS  code 
in  a  two  dimensional  setting.  The  interface  between  hard  (rigid,  two  dimensional) 
inclusions  embedded  in  an  elastomer  characterized  by  a  thnee-term  (rate  insensitive)  Ogden 
model,  is  represented  by  a  cohesive-zone  type  interfacial  model  to  follow  the  whole  process 
of  interfacial  dewetting  and  its  effect  on  the  global  (multiphase)  material  response  in  plane 
strain.  The  analysis  is  earned  out  through  a  mixed  finite  element  foiTnulation  for 
hyperelasticity,  incoiporating  interface  elements.  We  consider  the  effects  of  paiticle 
geometry  and  loading  conditions  on  the  process  of  interfacial  failure.  The  results  indicate 
that  the  distiibuted  failure  process  is  highly  unstable  and  depends  heavily  on  the  size, 
shape,  orientation  and  interactions  of  inclusions  as  well  as  the  global  loading  conditions. 
The  overall  material  behavior  of  the  model  agrees  qualitatively  with  experimental 
observation  this  instability  discovery  has  prompted  the  investigation  considered. 


Zhong,  X.A.  and  Knauss,  W.G; 

On  the  Stability  of  Phase  Separation  in  a  Finite  Solid  with  Interfaces*; 
Accepted  for  publication  in  Mechanics  of  Composite  Materials  and  Structures. 

Abstract 

The  stability  of  homogeneous  phase  separation  in  finite  solids  containing  planar,  cylindrical 
or  spherical  interfaces  is  investigated  analytically.  Explicit  stability  conditions  are  deduced 
for  each  interface  geometiy.  It  is  shown  how  the  interaction  of  load  (force  or  displacement) 
material  properties  of  die  phases  and  interface  properties  jointly  detennine  the  stability  of 
the  interface  separation  process. 


Zhong,  X.A.  and  Knauss,  W.G; 

Effects  of  Particle  Interaction  and  Size  Variation  on  Damage  Evolution  in 
Filled  Elastomers* 

Submitted  to  "Mechanics  of  Materials". 

Abstract 

A  micromechanical  analysis  of  damage  evolution  (interfacial  debonding)  in  particulate  filled 
elastomers  addresses  the  effect  of  the  interactions  between  particles  and  of  the  variation  in 
filler  size.  The  composite  is  treated  as  an  assembly  of  two  constituents  in  a  finite  element 
model.  The  interaction  between  paiticles  controls  the  damage  evolution:  (1)  For  high 
volume  fractions,  a  relatively  small  change  in  particle  size  has  a  suiprisingly  large  effect  on 
the  local  material  response,  (2)  for  lai-ge  differences  in  paiticle  sizes  {e.g  bimodal 
distribution),  damage  occurs  at  interfaces  between  large  particles  and  the  matrix,  with 
limited  damage  occuning  at  the  small  pai  ticles.  While  these  effects  of  paiticle  interaction 
and  size  variation  are  smoothed  out  in  a  lai'ge  ensemble  of  particles,  it  is  foreseeable  that 
they  are  an  important  factor  in  a  failure  process  such  as  macroscopic  crack  propagation, 
which  spans  scales  considerably  larger  than  the  maximum  particle  size.  Specifically,  one 
expects  thus  that  in  the  vicinity  of  a  macroscopic  crack  the  large  paiticles  become  the  sites 
for  small  cracks  which  coalesce  to  larger  ones  and  join  up  with  the  macro  crack,  while  the 
small  paiticles  operate  primaiily  so  as  to  locally  stiffen  the  matrix  without  incuning 
significant  damage  in  their  vicinity. 
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PARTICULATE  COMPOSITE 
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Graduate  Aeronautical  Laboratories,  California  Institute  of  Technology,  Pasadena, 

CA91125,U.S.A. 


ABSTRACT 

A  full  field  method  for  visualizing  strain  fields  around  a  crack  tip  under  large 
strains  is  developed.  Digital  image  correlation  is  used  to  compute  strains  and 
displacements  incrementally  between  consecutive  images  in  a  process  of  large 
deformations.  Values  of  strain  and  displacements  for  these  consecutive  deformations  are 
added  such  that  convergence  of  the  DIG  algorithm  is  assured.  The  method  is  used  to 
investigate  the  strain  distribution  in  a  globally  homogeneously  strained  particulate 
composite  (solid  propellant)  as  well  as  in  a  zone  close  to  (~2  mm)  the  crack  tip  in  such  a 
material  by  using  a  microscope.  It  is  found  that  maximal  strain  variations  deviate  by  as 
much  as  a  factor  of  three  from  the  average  strain;  additionally,  observations  on  the 
interaction  of  strain  inhomogeneities  with  the  tip  of  a  macroscopic  crack  are  discussed. 

Keywords:  A.  fracture,  A.  strain  localization,  B.  inhomogeneous  material,  B.  viscoelastic 
material,  C.  digital  image  correlation. 


1.  INTRODUCTION 


Particulate  composites  are  widely  used  in  engineering.  In  the  automotive  industry, 
for  example,  carbon  black  filled  rubbers  are  used  in  tires.  Many  injection  molded 
materials  are  filled  with  small  particles,  while  other  rigid  polymers  are  toughened  through 
the  addition  of  rubber  particles.  Solid  propellant  rocket  fuels  are  physical  mixtures  of 
mostly  ammonium  perclorate  and  aluminum  powder,  often  in  a  multimodal  size 
distribution,  bonded  together  by  a  rubber  phase  called  the  matrix.  Failure  in  all  these 
materials  is  heavily  dependent  upon  the  interaction  between  the  particles  and  the  matrix, 
specifically  on  the  separation  of  particles  and  binder.  Failure  also  depends  on  the  volume 
ratio  of  particles  to  matrix,  which  is  typically  close  to  75%  in  solid  propellant  materials, 
but  only  in  the  range  of  15  to  40%  in  stmctural  polymers.  In  the  sequel  we  examine  the 
failure  progression  in  a  solid  propellant  (Triokol  TPH  1011).  Application  of  continuum 
mechanics  to  the  stress/strain  analysis  of  structures  made  of  these  types  of  materials 
typically  invoke  macroscopically  homogeneous  material  performance,  even  though 
deformations  are  an5^hing  but  homogeneous  at  the  size  scale  of  the  particles.  We  shall 
see  that  inhomogeneous  deformations  occur  at  a  size  scale  that  is  significantly  larger  than 
the  largest  particle,  and  that  the  failure  process  is  directly  related  to  the  micro-structural 
deformations  associated  with  these  inhomogeneities. 

Measuring  large  strains  over  small  domains  of  tens  to  hundreds  of  microns  is  not 
a  trivial  matter.  Imprinted  grids  tend  to  serve  well  at  a  size  scale  just  above  what  is 
required  here.  Determining  the  micromechanical  deformation  with  the  aid  of  optical 
microscopy,  e.g.  at  the  tip  of  a  macroscopic  crack,  implies  the  need  to  extend  the 
presently  available  tools  of  strain  measurements.  In  principle  the  digital  image  correlation 
method  [Sutton  et  al,  1983,  1985,  1986  and  1989],  [Vendroux  /  Knauss,  1994]  is  ideal 
for  this  purpose  except  that  it  is  not  suitable  if  the  deformations  are  so  large  that 
convergence  of  the  correlation  algorithm  is  no  longer  guaranteed:  We  shall  see  that 
deformations  involving  strains  much  in  excess  of  10%  cause  convergence  failure  of  the 
Die  algorithm.  On  the  other  hand,  strains  on  the  order  of  50%  to  100%  are  typical  for 
crack  propagation  problems  in  the  materials  of  interest  here.  Accordingly  we  develop  an 
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incremental  application  of  the  DIG  method  that  is  capable  of  analyzing  large  deformation 
histories.  This  development  is  first  addressed  in  section  2,  and  followed,  in  section  3,  by 
a  discussion  of  the  experimental  setup  for  two  sets  of  experiments:  The  first  experiment 
addresses,  in  section  4,  locally  inhomogeneous  deformations  in  a  globally  homogeneous 
deformation  field,  and  the  second  examines  the  deformation  field  around  the  tips  of  a 
(slowly  growing)  crack  in  section  5,  with  particular  interest  centering  on  the 
inhomogeneity  of  the  material  response  in  the  immediate  crack  tip  region.  The  paper  is 
summarized  with  concluding  remarks  in  section  6. 

2.  DIGITAL  IMAGE  CORRELATION 

Developed  by  Sutton  and  his  colleagues  [1983,  1985,  1986  and  1989]  and 
improved  by  Vendroux  and  Knauss  [1994],  the  digital  image  correlation  (DIG)  program 
is  used  to  measure  the  displacement  field  and  its  gradients  from  images  of  an  undeformed 
and  deformed  body.  These  are  gray  level  images  consisting  of  a  grid  of  pixels,  (typically 
640  by  480)  with  eight-bit  gray  levels  (0  to  255  levels).  In  the  sequel  we  discuss 
problems  arising  with  large  deformations  and  a  remedy  to  the  situation  through  a  step¬ 
wise  method  we  call  Large  Deformation  Digital  Image  Gorrelation  (LD-DIG). 


2. 1.  Effect  of  Strain  Level  in  Code  Convergence 

The  problem  in  applying  the  DIG  program  to  compute  strain  distributions  in  a 
large  deformation  process  is,  essentially,  the  failure  of  the  DIG  algorithm  to  converge 
from  an  initial  solution  estimate.  The  reasons  for  non-convergence  may  be  diverse.  The 
two  major  ones  are  changes  of  shadows  from  a  fixed  light  source  coupled  with  large 
motions  of  the  surface,  and,  in  the  present  case,  possibly  an  inhomogeneous  evolution  of 
the  deformation  images.  If  one  considers  that  the  (rate  of)  convergence  depends  on  the 
closeness  of  an  initial  estimate  for  the  result  (see  the  study  on  the  radius  of  convergence 
by  Vendroux  and  Knauss  (1997)  it  becomes  reasonable  that  failure  can  occur  at  even 
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moderate  strains.  Clearly,  more  definitive  rules  for  code  failure  or  success  depend  on  the 
specific  experimental  conditions. 

To  assess  the  effect  of  the  strain  level  on  the  convergence  of  the  DIC  code, 
consider  a  test  on  a  homogeneous  silicone  rubber  sheet  stretched  uniaxially,  for  which  the 
resultant  undeformed  and  deformed  images  associated  with  stretches  from  0%  to  40%  are 
compared  with  the  aid  of  DIC.  For  each  deformation  the  strains  and  displacements  were 
computed  at  300  points.  The  fraction  of  points  at  which  the  numerics  for  the  correlation 
optimization  converged  is  presented  in  Figure  1  as  a  function  of  the  Lagrangian  strain. 


Figure  1.  Successful  convergence  of  the  DIC  algorithm  as  a  function  of  (Lagrangian) 

strain  level. 


It  is  apparent  that  for  deformations  in  excess  of  10%  a  pronounced  decrease 
occurs  in  the  number  of  points  with  successful  convergence.  For  the  purpose  of  studying 
cracked  solid  propellants,  where  typically  strains  in  excess  of  30%  need  to  be  measured, 
the  applicability  of  the  standard  DIC  method  is  thus  seriously  compromised  and  a  new  or 
extended  analysis  tool  is  required. 
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2. 1.  The  Large  Deformation  Digital  Image  Correlation  Method 


To  illustrate  a  proposed  Large  Deformation  Digital  Image  Correlation  (LD-DIC), 
consider  a  sequential  deformation  process  on  a  body:  Initially  undeformed,  the  body 
undergoes  a  continuous  deformation,  called  the  “global  deformation”.  Consider  three 
configurations  of  the  body  at  three  different  instants  during  this  global  deformation.  The 
first  configuration  describes  the  undeformed  state  of  the  body  and  the  second  a  deformed 
state  under  a  set  of  changed  surface  (and  body)  forces.  Call  this  first  segment  of  the 
global  deformation  “deformation  A”.  Next,  an  increment  in  the  surface  and  body  forces 
deforms  the  body  further,  this  next  incremental  deformation  being  designated  by  “B”. 
The  state  of  the  body  after  deformation  B  is  represented  by  the  third  configuration.  Each 
configuration  corresponds  to  an  experimentally  determined  and  temporally  ordered  set  of 
images,  say,  the  and  We  select  load  parameters  such  that,  by  assumption,  the 

Die  code  can  converge  successfully  to  the  proper  increments  of  displacements  and 
displacement  gradients  for  the  two  separate  deformations.  However,  the  strains  between 
configurations  1  and  3  (globiil  deformation)  are  presumed  larger  than  those  that  cause 
convergence  failure.  The  LD-DIC  method  computes  the  deformation  for  the  third  state 
from  the  deformations  in  the  two  separate  steps.  To  this  end  compute  the  global 
deformation  gradient  tensor  Fgiobai  as  the  product  of  the  individual  deformation  gradients 

Fglobal=rBFA  (1) 

The  Die  program  determines  the  displacements  of  deformation  A  and  their 
gradients  for  a  discrete  set  of  points  Gi  defined  on  a  rectangular  grid  with  respect  to  the 
first  (reference)  configuration.  For  demonstration  purposes,  these  points  are  represented 
in  configuration  1  in  Figure  2  as  a  coarse  rectangular  grid.  At  the  end  of  deformation  A, 

material  particles  of  the  body  at  the  grid  points  have  moved  to  the  points  Gi  in 
configuration  2  as  signified  by  the  non-orthogonal  grid.  By  comparing  the  position  of  the 
points  Gi  and  G  i  in  configurations  1  and  2  the  DIC  program  yields  the  displacements  \ 
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and  v'^  i,  and  the  associated  gradients  Ux^,  Vy\  Uy^  and  Vx^  where  the  index  “i”  signifies 
the  individual  initial  positions  of  the  Gi. 


Configuration  1  Configuration  2  Configuration  3 


Figure  2.  The  interpolation  process. 

The  displacements  and  displacements  gradients  are  computed  in  a  Lagrangian 
reference  frame  the  coordinates  of  which  are  X  =  {Xn  ,  X^}.  Upon  denoting  the 
displacement  components  as  the  motions  of  the  material  points  are  then  represented 
by 

G,(Xi)=  G,(Xi)+Ui^(Xi)  .  (2) 

During  the  second  deformation  (B),  the  displacements  and  their  gradients  are 
computed  for  a  different  set  of  material  particles  which  are  located  in  configuration  2  at 
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the  grid  points  Kj  of  an  orthogonal  grid,  shown  in  Figure  2  as  the  fine  grid.  In 

deformation  B  these  points  are  mapped  onto  the  set  Kj  in  configuration  3.  Kj  is  again 
represented  in  the  Lagrangian  reference  frame  Yi  =  {Yn  ,  Y2i}  of  configuration  2  as 


|.(Yi)=  ^/Yi)+u/(Yi)  (3) 

To  obtain  the  global  deformation  (A+B)  it  is  necessary  to  assure  that  the  material 
points  of  the  first  deformation  are  the  same  as  for  the  second  one.  It  is  thus  necessary  to 
interpolate  the  results  from  the  DIG  for  deformation  B  onto  the  locations  of  the  material 
particles  or  grid  points  of  the  first  (reference)  configuration  to  obtain  the  displacements 
and  their  gradients  u^i,  v®i,  u^i,  Vy®i,  Uy®i  and  Vx®i  relative  to  configuration  1.  The 
interpolation  at  the  coordinates  G,.  is  achieved  by  constracting  a  set  of  piecewise 
continuous  surfaces  from  bilinear  patches  (plane  surfaces).  These  are  determined  from 
the  coordinates  Kj  of  configuration  2  so  that  the  displacements  and  displacement 

gradients  of  the  initial  set  of  material  points  are  now  known. 

Using  the  subscript  “gl”  to  denote  global  variables,  one  finds  the  global 
displacement  components  as 

Ug,  +u® 

Vg,=v^+v®.  (4) 


and  their  gradients  are,  upon  invoking  the  tensorial  relation  (1), 


=«/ 

[v,]g,=v/ 

[vj,,  =v/ 


.  B  ,  A  B  ,  BA 


.  B  ,  A  B  ,  BA 

+  V  +V  V  +V  U 

*  y  y  y  X  y 


.  B  .  A  B  ,  BA 

+  My  +Mj,  +My  Vy 


.  B  ,  A  B  ,  BA 
+  V,  +U^  V,  +V^  V, 


(5) 
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From  which  the  (global)  two-dimensional  Lagrangian  strain  tensor  derives,  by 
definition,  as 


2Exx  —  2[Mx]gi  +  { [wjgi  +  [vjgi  } 

2Eyy  =  2[Vy]gl  +  {  [Wylgl  +  [Vy]g|  } 

2Exy  =  {  [My]g|  +  [vjgl  }+  {[Mjgl  [My]gl  +  [vjgl  [Vy]gl  }  (6) 

Here  the  displacement  gradients  of  the  displacement  “w”  normal  to  the  surface 
have  been  neglected.  This  is  permissible  (even)  for  this  situation  of  “plane  stress”  since 
under  these  large  deformations  the  propellant  material  voids  considerably  so  that  the 
effective  Poisson  ratio  is  rather  small,  thus  giving  rise  to  only  small  gradients.  For  a 
global  deformation  requiring  more  than  two  steps  the  same  is  applied  consecutively  such 
that  he  second-to-the  last  increment  is  treated  like  the  first  one  in  the  two-step  example 
outlined  here. 

When  more  than  two  steps  are  needed  for  field  evaluations  there  arises  a  loss  of 
information  at  the  boundary.  This  is  the  result  of  interpolating  information  near  and 
internal  to  the  boundary.  As  more  and  more  steps  are  required  the  information  near  the 
boundary  becomes  increasingly  corrupted  due  to  the  interpolation  process  and  thus 
information  is  lost  there.  This  loss  is  evident  in  the  field  images  presented  later  on  as 
apparent  white-out  regions  along  a  boundary,  including  the  flanks  of  a  crack. 
Minimization  of  this  feature  still  requires  future  attention. 

2. 3.  Verification  of  the  Scheme  for  Addition  of  Fields 

In  order  to  check  the  efficacy  of  the  proposed  multi-step  scheme  we  examine 
experimentally  a  specimen  of  homogeneous  silicone  rubber  without  a  crack  and  coated 
with  microscopic  speckles,  stretched  sequentially  and  uniaxially  to  a  maximum 
(Lagrangian)  strain  of  0.70  in  a  sequence  of  12  deformation  steps  of  3-4%  strain  each  (13 


8 


images).  These  strains  were  recorded  (optically)  with  the  aid  of  a  microscope  by  keeping 
track  of  special  markers  (=  prescribed  strain).  Also,  by  using  the  information  generated 
by  the  DIC  program  for  every  sub-deformation  and  the  Large  Deformation  DIG  method, 
the  strains  corresponding  to  images  1  and  2,  images  1  and  3,  etc.,  were  computed,  up  to 
the  deformation  of  image  15  relative  to  image  1. 


Figure  3.  Comparison  of  strains  determined  optically  and  by  the  LD-DIG  method;  the 

solid  line  is  the  ideal  relation. 

Figure  3  shows  the  result  as  a  comparison  of  the  (optically)  prescribed  strain  and 
the  Large  Deformation  DIC  computed  (Lagrangian)  strains.  The  maximum  deviation 
occurs  at  a  strain  of  40%,  amounting  to  only  a  1%  difference  between  the  strain 
determined  through  the  microscope  and  by  the  large  deformation  DIC  method,  a 
precision  that  is  very  acceptable  for  experimental  mechanics  investigations. 
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3.  EXPERIMENTAL  SETUP  FOR  DEFORMATION  STUDIES 


The  equipment  used  for  the  experimental  work  on  the  particulate  composite  (solid 
propellant)  has  been  described  elsewhere  in  detail  by  Gonzalez,  (1997).  Includes  a 
straining  stage  driven  by  a  stepping  motor  through  a  flexible  cable.  The  straining  stage  is, 
in  turn,  mounted  on  a  positioning  stage,  for  which  a  joy-stick  controller  allows  the 
positioning  of  the  straining  stage  under  the  objective  of  a  Nikon  microscope 
(Measurescope  MM-22),  a  CCD  camera  and  a  personal  computer  with  a  frame  grabber 
unit. 


FRAME 

GRABBER 


Figure  4.  Schematic  of  the  Experimental  Setup 


The  images  from  the  experiment  are  processed  on  a  Sun  workstation.  A  schematic  of  the 
experimental  setup  is  shown  in  Figure  5. 
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4.  LOCALLY  INHOMOGENEOUS  DEFORMATIONS  IN  A  GLOBALLY 
HOMOGENEOUS  DEFORMATION  FIELD. 


Before  considering  the  problem  of  the  strain  distribution  around  the  tip  of  a  crack 
we  examine  the  inhomogeneity  of  strains  in  a  solid  propellant  sample  subjected  to 
uniform  deformations  at  the  specimen  boundary  (globally  homogeneous  deformation).  To 
this  end  a  sheet  sample  of  propellant  [5  cm  x  2  cm  x  0.3  cm]  is  strained  uniformly  in  the 
direction  parallel  to  the  2  cm  dimension;  after  normalization  by  the  2  cm  dimension  that 
boundary  displacement  is  called  the  “applied  strain”.  In  Figure  5  we  show  a  2.5  mm  x 
2.55  mm  field  from  the  center  region  of  the  specimen  as  resolved  through  a  microscope 
and  deduced  with  the  aid  of  the  LD-DIC  algorithm.  The  false  color  scheme  clearly 
identifies  the  inhomogeneous  character  of  the  strain  field,  and  from  this  map  it  is  quite 
clear  that  the  variations  in  the  strain  values  are  not  only  very  significant,  but  that  the 
“material  properties”  vary  to  a  like  degree  in  these  regions. 

It  is  of  interest  to  dwell  briefly  on  the  scale  of  the  inhomogeneous  regions.  While 
the  latter  are  not  sharply  defined,  it  is,  nevertheless,  clear  that  these  domains  are 
measured  in  terms  of  millimeters  and  not  microns.  One  might  argue  that  a  variation  in 
properties  in  a  small  region  (a  hole  or  a  hard,  well-bonded  particle)  embedded  in  a 
homogeneous  field  renders  deviation  from  that  field  several  times  larger  than  the  defect 
itself  The  reference  problem  of  a  circular  hole  in  an  elastic  infinite  sheet  suggests  that  a 
noticeable  perturbation  in  its  vicinity  is  on  the  order  of  three  times  its  diameter.  While  it 
is  not  our  purpose  to  analyze  in  detail  here  the  precise  origin  of  every  inhomeogeneity,  be 
it  one  of  larger  strain  than  the  average  or  of  a  smaller  value,  it  appears  clear  that  these 
inhomogeneities  are  produced  either  by  clusters  of  particles  or  by  the  presence  or  absence 
of  individual  ones. 

To  quantify  the  inhomogeneity  of  the  strain  field  further,  consider  a  plot  of  the 
strain  along  the  line  in  Figure  5  as  shown  in  Figure  6.  The  magnitude  of  these  strains  vary 
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by  as  much  as  a  factor  of  three,  although  the  average  of  the  strain  in  Figure  6  represents 
closely  the  applied  global  strain  of  1.5  %. 

The  distribution  of  these  inhomogeneities  defines  a  macroscopic  size  scale  below 
which  the  assumption  of  homogeneous  material  properties  is  not  justified.  Stated 
alternately,  in  order  to  be  able  to  assign  homogeneous  properties  to  such  a  composite  it  is 
necessary  to  deal  with  a  size  scale  that  is  several  times  larger  than  the  spacing  between 
the  regions  of  inhomogeneity.  That  region  is,  however,  at  least  on  the  order  of  five  mm  or 
more. 


0.04 

0.035 

0.03 

0.025 

g  0.02 

0.015 

0.01 

0.005 


>  '  ■  I  ■  I  I  I  I  I  I  I  I  ■  1  I  I  I  1 

0.5  1  1.5  2  2.5 


Y(miin) 

Figure  6.  Eyy  Distribution  along  the  trace  x=  1.4mm  in  Figure  5,  parallel  to  the  tension 

axis. 
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We  shall  see  that  this  size  scale  is  important  in  the  analysis  of  the  failure/fracture 
process  considered  next,  since  in  this  context  this  size  limitation  emphasizes  the 
dependence  of  the  fracture  process  on  the  statistical  nature  of  the  medium. 

5.  THE  STRAIN  FIELD  NEAR  THE  TIP  OF  A  CRACK 

Having  examined  the  distribution  of  strains  in  a  globally  homogeneous 
deformation  field  we  turn  next  to  examining  the  deformations  in  the  close  vicinity  of  a 
crack  tip.  We  describe  first  the  experimental  set-up  and  then  proceed  to  the  analysis  of 
the  results. 

5. 1.  Experimental  Aspects 

A  cracked  specimen  of  solid  propellant  TPH  1011  is  deformed  globally  at  a 
constant  strain  rate  0.001  1/sec  in  the  direction  perpendicular  to  the  crack.  The  crack,  cut 
initially  with  a  razor  blade,  opens  commensurately.  Its  opening  process  is  monitored 
through  a  microscope  at  25x  power.  Five  digital  images  of  640  x  480  pixels,  representing 
4  mm  X  5  mm  of  the  specimen  surface  were  acquired  every  10  seconds,  corresponding  to 
global  (Lagrangian)  strains  of  Eyy  =0%,  1%,  2%,  3%  and  4%. 


Figure  7.  Specimen  geometry  and  straining  of  a  specimen. 
The  dotted  lines  represent  the  deformed  geometry. 
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Figure  5:  Inhomogeneous  strain  distribution  in  an  integral  specimen  “homogeneously” 

deformed  to  1.5%  Global  Strain’. 


X  (mm) 


Figure  8:  Maximum  Principal  Strain  Distribution  for  2%  Global  Strain. 


'  Note  that  the  red  domains  surrounding  the  (green)  high  strain  regions  are  the  result  of  clustering  of  the  red 
contour  lines,  not  strain  concentrations  in  themselves. 
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Although  each  of  the  images,  excepting  the  first,  is  associated  with  different 
deformations  we  select  here  only  two  loads  or  deformation  levels  for  illustrative 
purposes.  Application  of  the  (multi-step)  LD-DIC  code  renders  deformation  maps  as 
shown  in  figures  8  and  9. 

Recall  that  because  the  multi-step  method  looses  information  at  the  boundary  in 
each  step  the  crack  appears  wide  open  and  as  having  a  very  blunt  tip.  This  appearance  is 
the  consequence  of  the  multi-step  method,  so  that  Figures  8  and  9  do  not  represent  the 
crack  opening  shape  correctly;  instead  the  shape  has  been  sketched  in  as  outlines. 

5. 2.  Results 

We  demonstrate  results  for  the  global  strain  of  2  and  4%.  While  deformations  at 
higher  load  levels  can  be  obtained,  their  interpretation  is  more  troublesome  since  it 
involves  (more)  motion  of  the  crack  tip.  As  an  example  we  present  (false  color)  plots  of 
the  maximum  principal  strain  at  2  %  global  deformation  in  Figure  8  and  both  the 
maximum  and  minimal  principal  strains  for  4  %  in  Figures  9a  and  9b.  The  two  strain 
levels  are  presented  primarily  to  afford  a  comparison  for  two  progressively  larger  strains. 
By  comparing  Figures  8  and  9a  one  can  readily  see  how  the  inhomogeneities  develop 
early  and  essentially  grow  in  intensity  with  the  global  strain.  In  passing  from  2  to  4  %  of 
global  strain,  a  small  amount  of  crack  growth  has  taken  place  as  is  evident  from 
comparing  Figures  8  and  9a. 

As  before,  the  amplitude  of  the  (Lagrangian)  strain  is  represented  by  colors  and 
contours  and  the  small  lines  represent  the  orientation  of  the  corresponding  principal  axes. 
Note  that  although  the  maximum  principal  axes  should  be  basically  parallel  to  the  applied 
global  displacement(s),  at  least  in  the  region  ahead  of  the  crack  tip,  there  are  numerous 
locations  where  marked  differences  occur  from  this  orientation.  These  differences  are 
associated  with  the  material  strain  inhomogeneities  discussed  in  section  4.  The  same 
observation  applies  to  the  map  of  the  minimum  principal  strain  in  Figure  9b  for  which  the 
orientations 
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Figure  9.  Principal  Strain  Distribution  for  a  Far  Field  Strain  of  4% 

(a)  Maximum  Principal  Strain  (upper) 

(b)  Minimum  Principal  Strain  (lower) 
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(small  line  segments)  are  generally  orthogonal  to  those  for  the  maximum  principal  strain. 
Also,  compare  the  scale  associated  with  the  inhomogeneities  and  the  distance  of  the  high 
strain  region  relative  to  the  crack  tip.  It  is  clear  that  the  strain  inhomogeneities  cannot  be 
separated  from  or  submerged  in  the  strain  field  of  the  crack  tip:  The  crack  tip  strains  are 
intimately  connected  with  the  scale  of  the  material  inhomogeneity.  We  note  in  particular 
that  at  this  size  of  the  viewing  field  the  lobes  that  are  part  of  the  typical  crack  tip  strain 
field  for  isotropic  and  homogeneous  material  are  absent  in  the  domain  of  observation  or 
at  most  apparent  in  vestigial  form;  they  are  present  only  outside  of  the  2.5  mm  x  2.5  mm 
observation  area. 

One  feature  of  interest  and  identifiable  in  figures  8  and  9  is  the  localization  of 
high  strain  within  a  roughly  circular  area  of  0.5  mm  radius  around  the  crack  tip.  Figure  8, 
which  shows  the  maximum  principal  strains  and  directions  for  2%  far  field  strain,  shows 
deformation  localization  around  the  crack  tip  in  two  locations.  One,  partially  visible,  is 
centered  on  the  current  crack  tip  (at  x=0.7mm,  y=l.lmm)  and  about  150  microns  in 
radius.  The  other  concentration  is  centered  on  the  position  (0.7mm,  0.8mm),  where  it 
reaches  a  maximum  principal  strain  of  about  11%.  However,  that  location  does  not 
become  a  part  of  the  crack  propagation  process  as  is  evident  by  its  persistence  in  Figure 
9a.  Adjacent  to  these  strain  concentrations  is  a  domain  of  about  5%  local  strain  which 
includes  the  position  (1.1mm,  1mm)  where  the  orientation  of  the  maximum  principal 
strains  is  nearly  aligned  with  the  crack  rather  than  being  normal  to  it^.  At  this  position  a 
void  develops  under  subsequently  increased  deformation.  This  void  is  visible  in  Figure 
10,  which  represents  the  image  of  the  specimen  for  a  5%  global  deformation. 

We  have  dwelled  on  these  details  to  some  extent  in  order  to  demonstrate  that  the 
observation  of  the  surface  of  the  failing  material  is  not  always  indicative  of  how  and 
where  the  crack  is  likely  to  propagate.  Because  one  identifies  strain  inhomogeneities  on 
the  surface,  it  stands  to  reason  that  similar  distributions  exist  through  the  thickness  of  the 

^  Similar  situations  prevail  in  Figure  7  at  locations  (1.6  mm,  0.9  mm),  (2.1  mm,  1.1  mm)  and  (2.2  mm,  1.6 
mm).  In  these  regions  the  strains  remain  small,  indicting  that  they  are  associated  with  rigid  domains  inside 
the  specimen  and  under  the  surface. 
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specimen.  Thus,  the  observations  offered  here  point  to  a  truly  three-dimensional  process. 
However,  it  is  quite  clear  that  the  domain  in  which  the  failure  process  is  prominently 
operative  is  confined  to  a  domain  on  the  order  of  half,  but  not  more  than  one  millimeter. 
This  observation  agrees  with  the  results  of  Liu  [Liu,  C.T.,  1997]  who  observed  a  similarly 
confinement  of  the  process  zone  to  a  very  small  region.  So  much  is  clear  from  these 
studies,  only  portions  of  which  are  reported  here,  namely  that  crack  propagation  occurs 
by  opening  up  voids  which  are  typically  high  strain  regions,  distributed  statistically 
throughout  the  material,  and  through  joining  of  these  voids  with  the  main  crack. 


Figure  10.  Micrograph  of  the  crack  tip  region  at  a  global  strain  of  5%.  The  “void”  at  the 
crack  tip  was  a  region  of  “small”  (~5%)  surface  strain  in  Figure  8. 
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6.  CONCLUDING  REMARKS 


The  method  of  digital  image  correlation  has  been  extended  to  large  deformations 
by  dividing  the  strain  range  into  intervals  within  each  of  which  DIG  converges.  The 
sequential  extension  renders  satisfactory  results  with  deviations  not  exceeding  1%  from 
the  prescribed  strain  field.  This  method  has  allowed  resolution  of  highly  inhomogeneous 
deformations  embedded  in  a  globally  homogeneous  deformation  field.  It  is  then 
demonstrated  that  this  method  can  be  applied  to  the  analysis  of  strains  around  the  tip  of  a 
crack  in  a  particulate  composite  (solid  propellant  rocket  fuel). 

The,  perhaps,  most  striking  result  of  that  investigation  is  the  surprisingly  large, 
inhomogeneous  variations  in  the  local  strain  field  in  a  globally  homogeneously  deformed 
solid  propellant  material.  It  is  apparent  that  the  inherent  heterogeneity  of  the  material 
plays  a  key  role  for  distribution  of  strains  around  the  crack  tip  and  its  propagation.  A 
second  important  observation  is  the  fact  that  most  of  the  deformation  related  to  the  crack 
propagation  process  localizes  in  a  region  around  the  crack  tip  of  only  about  0.5mm 
radius.  Although  there  exist  strain  concentration  domains  outside  of  this  small  region  on 
or  close  to  the  specimen  surface,  these  deformations  are  insufficiently  high  to  cause  void 
formation  of  the  strain  concentrations  so  that  they  do  not  become  sources  of  coalescence 
with  the  macroscopic  crack.  There  is  evidence,  however,  that  the  three  dimensional 
distribution  of  strain  inhomogeneities  through  the  thickness  of  the  specimen  plays  a  role 
that  mimics  their  in-plane  distribution.  There  is  no  reason  to  suppose  otherwise  for  this 
kind  of  particulate  solid.  Certainly  some  of  the  features  visible  on  the  surface  are  the 
consequence  of  single  particles  and  particle  agglomerations  buried  just  beneath  the 
surface. 
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Analysis  of  Interfacial  Failure  in 
Particle-Filled  Elastomers 

The  evolution  of  microdamage  (interfacial  dewetting)  in  highly  filled  elastomers 
under  consideration  of  high  deformation  gradients  is  examined.  The  interface  between 
hard  (rigid,  two-dimensional)  inclusions  embedded  in  an  elastomer  characterized  by 
a  three-term  Ogden  ( rate  insensitive )  model,  and  the  elastomer  matrix  is  represented 
by  a  cohesive-zone  type  interfacial  model  to  follow  the  whole  process  of  interfacial 
dewetting  and  its  effect  on  the  global  (multiphase)  material  response  in  a  plane 
strain  setting.  The  analysis  is  carried  out  through  a  mixed  finite  element  formulation 
for  hyper  elasticity,  incorporating  interface  elements.  We  consider  the  effects  of  parti¬ 
cle  geometry  and  loading  conditions  on  the  process  of  interfacial  failure.  The  results 
indicate  that  the  distributed  failure  process  is  highly  unstable  and  depends  heavily 
on  the  size,  shape,  orientation  and  interactions  of  inclusions  as  well  as  the  global 
loading  conditions.  The  overall  material  behavior  of  the  model  agrees  qualitatively 
with  experimental  observation. 


1  Introduction 

Nonlinear  behavior  is  widely  observed  in  particulate  compos¬ 
ites,  such  as  solid  propellants,  tires,  toughened  plastics,  even 
when  global  deformations  are  relatively  small.  Global  non-lin¬ 
earity  in  particulate  composites  is  often  caused  by  cavitation  or 
tear  (appearance  of  voids  in  the  matrix)  and/or  by  interfacial 
particle  debonding  or  “dewetting”;  see  Schippel  (1920),  Smith 
(1959),  Fanis  (1964),  Oberth  (1967),  Knauss  and  Mueller 
(1979),  and  Gent  and  Park  (1984).  Dewetting  and  interfacial 
void  generation  occurs  when  the  bond  at  an  interface  is  rela¬ 
tively  weak;  otherwise  cavitation  can  occur.  How  to  represent 
this  damaged-material  behavior  in  continuum  terms  for  finite 
element  analyses  has  been  a  long  standing  question.  Constitutive 
models  containing  damage  have  been  proposed,  for  example, 
by  Farris  and  Schapery  (1973),  Schapery  (1986),  Govindjee 
and  Simo  ( 1992),  Vratsanos  and  Farris  (1993),  and  Ravichan- 
dran  and  Liu  ( 1995).  However,  these  models  (except  Govindjee 
and  Simo,  1992)  are  typically  restricted  to  infinitesimal  elastic/ 
viscoelastic  matrix  behavior,  and  therefore  their  applications  in 
a  high  deformation  gradient  area,  such  as  a  crack  tip,  are,  at  best, 
questionable.  More  importantly,  one  must  critically  examine  the 
implied  proposition  of  first  representing  a  truly  discontinuous 
process  by  a  homogenized  continuum  formulation  and  then  turn 
around  and  formulate  a  new  local  failure  criterion  that  addresses 
the  crack  propagation  process  by  means  of  the  homogenized 
material  description  in  such  a  way  that  the  true  physical  process 
of  the  discrete  micro  crack  growth  and  coalescence  is  repro¬ 
duced  faithfully  near  the  crack  tip. 

By  contrast,  it  appears  more  reasonable  to  characterize  the 
material  locally  near  a  crack  tip  in  a  discrete  way  by  modeling 
a  select  region  of  the  high  strain  gradient  domain  through  dis¬ 
cretely  failing  elements  embedded  in  another  and  larger  region 
of  nonlinear  and/or  linear  material  response  as  illustrated  in 
Fig.  1.  Such  a  hybrid  discrete-continuum  approach  appears  fea¬ 
sible  in  light  of  the  ever  more  rapidly  increasing  power  of 
computing  machines.  We  thus  pursue  here  initially  the  time- 
independent  (nonviscoelastic)  problem  of  the  dewetting  process 
allowing  for  large  deformation  in  the  matrix  material,  and  with 
the  intention  of  incorporating  the  results  into  a  larger  finite 
element  analysis  for  crack  growth  studies  in  particulate  compos¬ 
ites. 
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In  addition  to  interfacial  debonding,  we  consider  here  the 
formation  and  growth  of  voids,  as  well  as  their  interaction  in 
(high)  strain  gradients.  The  aim  is  to  apply  these  physical  mi¬ 
cro-phenomena  to  global  failure  in  these  types  of  composites: 
macroscopic  regions  of  stress  concentrations  are  to  be  repre¬ 
sented  by  an  assemblage  of  micro-mechanically  detailed 
“superelements,”  their  individual  force  transmissions  and  fail¬ 
ure  responses  governing  the  failure  procession  of  global  com¬ 
posites.  Thus  our  interest  is  not  limited  to  unit-cell  problems, 
but  aims  at  much  larger  scale  problems  in  order  to  address  the 
failure /fracture  process  in  particulate  composites.  Because  the 
domain  for  this  distinctly  discrete  process  of  failure  has  been 
shown  to  be  very  limited  in  extent  (Liu,  1991;  Gonzalez  and 
Knauss,  1997)  it  may  suffice  to  deal  with  a  rather  small  number 
of  “superelement,”  say  on  the  order  of  200  to  300. 

To  this  end  one  needs  to  characterize  the  debonding  of  parti¬ 
cles  from  the  matrix  material  with  the  key  issue  being  the  char¬ 
acterization  of  the  properties  of  a  matrix-particle  interface,  in¬ 
cluding  its  failure  response.  This  issue  has  been  addressed  re¬ 
peatedly,  for  example,  by  Ungsuwarungsri  and  BCnauss  ( 1987), 
Needleman  (1987, 1990),  and  Yeh  ( 1992).  Interfacial  constitu¬ 
tive  models  are  sometimes  called  cohesive-zone  type  models 
as  they  resemble  the  cohesive-zone  model  (Barenblatt,  1962; 
Dugdale,  1960)  in  fracture  mechanics.  There  are  basically  two 
approaches  to  characterize  an  interface  at  the  continuum  level. 
One  utilizes  a  nonlinear  spring  model  to  mimic  the  interfacial 
failure  process;  the  other  employs  a  potential  for  formulating 
an  inteifacial  “constitutive  relation.”  The  common  feature  is 
that  they  relate  interfacial  tractions  to  relative  displacements  at 
interfaces.  Here  we  employ  a  simple  nonlinear  foundation 
model  which  is  similar  to  that  of  Ungsuwarungsri  and  Knauss 
( 1987),  the  parameters  of  which  can  be  ultimately  determined 
by  fitting  numerical  solution  to  corresponding  experimental  re¬ 
sults.  We  implement  the  proposed  interfacial  constitutive  rela¬ 
tion  through  interface  elements  which  have  been  widely  used 
in  computational  contact  problems,  rock  joints,  concrete  me¬ 
chanics;  see  for  example.  Beer  ( 1985),  Zubelewicz  and  Bazant 
(1987),  and  Schellekens  and  De  Borst  (1993). 

Though  an  elastomer  is  a  time-dependent  viscoelastic  mate¬ 
rial,  we  consider  the  elastomer,  for  simplicity  reasons,  to  be 
characterized  by  a  hyperelastic  Ogden  material  (Ogden,  1972). 
Because  the  elastomer  matrix  is  soft  compared  to  the  particles, 
we  consider  the  latter  to  be  rigid,*  In  order  to  gain  insight  into 


‘  This  choice  is  made  for  convenience.  Under  high  rate  loading,  the  binder 
may  be  stiff  by  compression  so  that  the  particles  will  also  deform.  That  case  is 
computationally  treatable  as  well,  but  requires  larger  computational  resources. 
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Fig.  1  A  discretely  failing  region  embedded  in  a  linear/nonlinear 
material 


this  approach  and  to  test  numerical  strategies,  a  two  dimensional 
version  of  the  problem  is  considered  first.  Furthermore  the  bond 
at  interfaces  is  assumed  to  be  “weak”  so  that  only  interfacial 
debonding  needs  to  be  considered  (no  cavitation  in  the  elasto¬ 
mer).  A  three-dimensional  problem  description  is  left  to  future 
investigations;  it  is  more  demanding  only  in  computational 
power. 

The  numerical  analysis  of  debonding  at  interfaces  has  been 
carried  out  by  several  authors,  for  example,  Lo  et  al.  (1994), 
Needleman  (1987,  1990),  Narasimhan  (1994),  Xu  and  Nee- 
dleman  (1994),  Yeh  (1992)  through  two  basic  approaches: 
One  is  the  so-called  nodal  release  technique  (Lo  et  al,  1994, 
Narasimhan,  1994),  the  other  employs  interface  elements.  The 
nodal  release  technique,  though  more  frequently  used,  permits 
only  a  single  node  to  be  released  at  a  time  by  reducing  the 
reaction  force  at  the  crack  tip  node  over  several  time  steps  and 
requires  continual  external  monitoring  and  interruption  of  the 
regular  finite  element  computation.  In  addition,  the  crack  propa¬ 
gation  path  or  debonding  path  has  to  be  prescribed  or  tracked 
node  by  node  which  makes  the  application  of  the  techniques  in 
cases  such  as  debonding  at  multiple  particles  rather  complicated. 

Interfacial  constitutive  relations  and  interface  element,  on  the 
other  hand,  can  treat  debonding  at  interfaces  automatically.  The 
present  approach  is  similar  to  that  of  Needleman  (1987).  The 
main  numerical  challenges  in  the  present  work  are  that  ( 1 )  the 
interfacial  failure  process  in  filled-elastomers  is  usually  less 
stable  than  that  in  other  types  of  particulate  composites  because 
of  the  relatively  soft  matrix  (Zhong  and  Knauss,  1997)  and  (2) 
The  local  instability  together  with  very  large  strains  occur  at  the 
interfaces  (tens  to  several  hundred  percent)  raises  convergence 
questions. 

In  the  sequel,  we  delineate  in  Section  2  the  interface  model 
along  with  the  associated  stability  issue.  The  computational 
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model  is  presented  in  Section  3  including  the  FEM  formulation 
for  hyperelasticity  and  the  interface  element.  Numerical  results 
are  demonstrated  in  Section  4  for  unit  cell  problems,^  to  evalu¬ 
ate  effects  of  inclusion  (particle)  size,  shape  and  orientation, 
as  well  as  of  boundary  loading  and  constraints  on  the  response 
of  a  unit  cell.  In  solid  propellants,  particle  sizes  vary  randomly 
in  a  range,  and  they  are  statistically  distributed  (Gonzalez  and 
Knauss,  1997).  To  account  for  this,  we  also  analyze  multi-cell 
problems,  considering  interactions  among  inclusions  of  identi¬ 
cal  or  different  sizes.  Additional  numerical  issues  are  discussed 
in  Section  5  ending  with  conclusions  in  Section  6. 

2  An  Interfacial  Constitutive  Model 

To  imitate  the  failure  behavior  of  cohesive  interface  material, 
the  interfacial  traction  is  assumed  to  increase  with  increasing 
interfacial  separation  to  a  critical  value,  after  which  the  traction 
decreases  with  increasing  separation  and  eventually  vanishes, 
corresponding  to  total  separation.  We  take  here  a  damage  me¬ 
chanics  point  of  view  by  assuming  that  the  interface  will  lose 
its  stiffness  as  interfacial  separation  reaches  this  critical  value 
(cf  Ungsuwarungsri  and  Knauss,  1987). 

2.1  The  Interfacial  Constitutive  Model.  Let  n  be  the 
normal  to  an  .interface,  and  s,- 1  the  corresponding  orthogonal 
tangential  directions;  E„ ,  and  E,  are  elastic  moduli  of  a  nonlin¬ 
ear  (dissipative)  spring  in  the  directions  (n,  s,  t).  The  spring 
response  is 

^  ES  0  0  “ 

D=  0  £?  0  ,  (1) 

0  0  £? 


as  long  as  the  normal  separation  Aw„  is  less  than  the  critical 
normal  separation  >^en  Aw„  >  Am^,  the  spring  moduli 
degrade  according  to 


(2) 

for  Am^  <  Am„  <  2Awc»  and  vanish  for  Am„  >  2Amc* 

The  nonlinear  interface  constitutive  relation  is  thus 

f  =  DAU(„,j,  (3) 

with  Am  the  relative  displacement  in  local  coordinates  (n,  s, 
t)  and  f  the  conjugate  interfacial  traction  defined  in  the  unde¬ 
formed  reference  configuration. 

Only  failure  due  to  normal  separations  is  accounted  for  in 
the  current  model.  Obviously  one  needs  to  consider  shear  failure 
at  an  interface  for  more  general  situations.  It  is  not  difficult  to 
include  shear  failure  into  the  above  interfacial  model;  see  for 
example  Tvergaard  and  Hutchinson  (1993).  However,  for  the 
type  of  problems  under  consideration,  the  shear  strength  of  the 
interface  is  much  larger  than  the  tensile  strength,  so  that  we  do 
not  consider  the  shear  mode  failure  in  the  present  work. 

2.2  Stability  of  an  Interfacial  Failure  Process  Suo  et 
al.  (1992)  investigated  the  stability  of  solids  containing  inter¬ 
faces  by  considering  two  semi-infinite  solids  bonded  along  a 

^We  use  this  tenn  frequently  to  refer  to  1 -inclusion  problems  (in  finite  do¬ 
mains)  for  notation  convenience,  even  though  the  analysis  does  not  always  impose 
periodic  boundary  conditions  that  belong  to  a  real  cell. 
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planar  interface  by  defining  the  instability  of  the  system  in  terms 
of  the  onset  of  nonunique  solutions  to  the  associated  boundary- 
value-problem.  Zhong  and  Knauss  ( 1997)  have  recently  studied 
the  stability  of  finite  solids,  with  planar  or  circular  interfaces 
under  homogeneous  decohesion  (i.e.,  the  whole  interface  de¬ 
bonds  uniformly  at  the  same  time),  deriving  stability  criteria  for 
a  general  form  of  the  traction-displacement  interfacial  relation 

/  =  /(Am).  (4) 

For  homogeneous  decohesion  at  a  bi-material  planar  interface 
the  stability  condition  is 

(5) 

where  corresponds  to  total  separation  at  the  interface.  The 
right-hand  side  of  inequality  (5)  represent  the  maximum  value 
of  -{dfld^u)  when  Am  is  in  (0,  Am,).  It  is  clear  that  the  size 
of  the  specimen  (L),  effective  elastic  modulus  =  E1E2I 
(El  -h  £2))^  where  Ei  and  E2  are  the  elastic  moduli  of  the  two 
materials)  and  the  unloading  part  of  the  interfacial  relation  (4) 
determine  the  stability  of  the  system.  A  similar  stability  condi¬ 
tion  is  obtained  for  a  circular  interface  (see  Zhong  and  Knauss, 
1997). 

An  interfacial  failure  process  associated  with  inhomogeneous 
decohesion  at  interfaces  is  intuitively  more  stable  than  that  asso¬ 
ciated  with  a  homogeneous  decohesion.  This  is  so  because  the 
surrounding,  bonded  regions  tend  to  stabilize  the  decohesion 
process.  Thus  the  stability  conditions  obtained  by  Zhong  and 
Knauss  (1997),  though  applied  only  to  homogeneous  decohe¬ 
sion  at  interfaces,  provide  conservative  estimates  for  the  onset 
of  instability  in  an  interfacial  failure  process.  It  is  also  clear 
from  the  above  stability  condition  (5),  that  the  more  compliant 
the  matrix  material,  the  less  stable  the  corresponding  composite 
is  under  interfacial  failure.  Thus  particulate-filled  elastomers 
tend  to  be  less  stable  than  particulate  metallic  composites  with 
respect  to  the  interfacial  failure  process. 

It  should  be  noted  here  that  the  instability  due  to  interfacial 
decohesion  is  local  because  of  the  sudden  drop  of  tractions  at 
the  interfaces.  This  local  instability  poses  a  problem  for  the 
numerical  simulation  in  a  quasi-static  setting  because  it  can 
(and  does)  lead  to  the  divergence  of  numerical  computations. 

3  Formulation  of  the  Computational  Model 

We  treat  the  matrix  as  a  hyperelastic  material  characterized 
by  a  3-tenn  Ogden  strain  energy  function  for  finite  deforma¬ 
tions.  For  reference  purpose,  we  summarize  the  hyperelasticity 
theory  and  its  hybrid  FEM  formulation  briefly. 

3.1  Hyperelasticity  and  Its  FEM  Formulation.  Let 
be  the  undeformed  reference  configuration  of  the  hyperelastic 
material (s),  the  current  deformed  configuration,  and  Vq  and 
Vi  the  corresponding  volumes.  XeQoi  y^H,.  Then 

y  =  X  +  u(X),  F  =  ^,  /=det(F).  (6) 

aX 

Here  u(X)  is  the  displacement  vector  at  material  point  X,  F 
the  corresponding  deformation  gradient  and  J  the  unit  volume 
change. 

Because  the  hyperelastic  material  is  (almost)  incompressible, 
the  deformation  gradient  F  is  multiplicatively  decomposed  into 
its  deviatoric  and  dilatational  parts  for  purpose  of  the  FEM 
formulation. 


The  strain  energy  function  for  a  nearly  incompressible  hyper¬ 
elastic  material  is  expressed  as  (Ogden,  1984), 

W(F)  =  W(F^n  +  (8) 

where  W  accounts  for  the  deviatoric  deformation  and  <{)  with 
<^(  1 )  =  0,  for  dilatational  deformation. 

Considering  frame  indifference,  and  further  assuming  that 
the  hyperelastic  material  is  isotropic,  we  have 

W(F)  =  u;(\i,\2,X3)  +  <^(7),  (9) 

where  X,  are  the  principal  deviatoric  stretches. 

Ogden  ( 1972)  assumes  that  u  has  the  form  (1972), 

w(\..  u,  \3)  =  I  -T  -  3)  (10) 

N 

SO  that  the  initial  shear  modulus  is  ^  =  2  fii  and  the  relation 
between  Cauchy  stress  and  the  principal  stretch  is 


If  the  material  is  fully  incompressible  (7  =  1 ),  we  have  (^(7) 
=  0,  and  the  corresponding  stress-stretch  relation  is 


Ti 


Pi 


(12) 


where  P,-  is  an  arbitrary  pressure  and  X,  are  the  principal 
stretches. 

For  almost  incompressible  material,  the  usual  displacement 
finite  element  formulation  can  behave  poorly  because  the  effec¬ 
tive  bulk  modulus  is  very  large  compared  to  its  shear  modulus, 
which  cause  the  stiffness  matrix  to  be  almost  singular  from  a 
numerical  point  of  view.  Consequently,  the  stress  calculated  at 
the  numerical  integration  points  show  large  oscillation  in  the 
pressure.  To  avoid  this,  a  hybrid  element  or  mixed  formulation 
has  been  proposed  (see  for  example,  Simo  and  Taylor,  1991). 
Although  extra  (pressure  or  compressibility)  variables  are  intro¬ 
duced  in  the  mixed  formulation,  the  effective  degree  of  freedom 
at  each  node  in  the  final  FEM  formulation  are  determined  by 
displacements  only,  because  the  extra  variables  are  eliminated 
at  the  element  level. 


3.2  Interface  Element.  The  interfacial  constitutive  rela¬ 
tion  is  implemented  by  constructing  an  interface  element  via 
the  principle  of  virtual  work  (Beer,  1985):  (1)  Set  up  local 
coordinates  (n,  s,  t)  to  describe  the  interface  geometrically. 
(2)  Find  the  relative  displacements  at  the  interface  in  local 
coordinates.  (3)  Use  the  principle  of  virtual  work  and  the  in¬ 
terfacial  constitutive  relation  to  obtain  an  equilibrium  equation 
for  the  interface  element  in  terms  of  global  nodal  coordinates. 
One  then  assembles  the  interface  element  stiffness  matrix  with 
remaining  element  stiffness  matrixes  to  form  the  global  stiffness 
matrix. 

Drawing  on  the  interfacial  constitutive  relation  and  the  princi¬ 
ple  of  virtual  work,  we  have 

f  B^fds  -  E%  (13) 

V  interface 

or 

f  B^DBads  =  F',  (14) 

V  interface 


with  I  the  identity  tensor.  It  is  easy  to  check  that  det  (F'^®'')  = 
1,  det  (F'"®^)  =  7  so  that  is  the  volume  preserving  part  of 
F,  and  F''°^  its  dilatation  part. 


with  the  force  vector  for  nodes  of  the  interface  element, 
a  is  the  element  displacement  vector  and  B  the  transformation 
matrix  between  nodal  displacements  and  displacement  in  the 
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local  coordinates.  The  residual  force  due  to  the  interface  ele¬ 
ment  is 

(15) 

interface 

and  the  tangent  stiffness  matrix  of  the  interface  element  is 


K,=  - 


m 

du  * 


(16) 


Here  n  refers  to  time  steps,  k  refers  to  the  sequence  of  iterations. 

For  the  interfacial  relation  (3),  the  tangent  stiffness  matrix 
for  interface  elements  is  nonsymmetric.  Thus  using  the  tangent 
stiffness  matrix  may  be  a  source  of  numerical  difficulty. 

We  use  direct  iteration,  which  is  expressed  for  the  interface 
elements  by 

f  B^D^„Bds[Aa] 

interface 

=  f  B^D!:,Bds{aU  (17) 

V  interface 

where  =  +  Aa,  n  is  the  time  step,  k  is  the  iteration 

index.  The  implementation  of  the  FEM  formulation  for  hyper¬ 
elasticity  is  available  in  ABAQUS  5.5.  We  incorporate  the  inter¬ 
face  element  (interfacial  constitutive  model)  into  ABAQUS  via 
a  user  subroutine  (user  element).  It  should  be  noted  that  the 
interface  element  has  zero  thickness  initially. 


4  Numerical  Results 

We  consider  plane  strain  deformation  of  the  “particulate” 
composites.  By  preliminary  choice  the  parameters  in  the  Ogden 
function  are  determined  from  fitting  experimental  data  of  Tre- 
loar  ( 1940^.  The  parameters  are  =  58.23,  aj  =  1.029,  ^2  = 
1.35  X  10  ^  a2  =  10.707,  ^3  =  0.246,  aj  =  -2.957.  They 
satisfy  the  material  stability  condition.  The  parameters  in  the 
interfacial  model  are  also  determined  quite  arbitrarily  for  now, 
except  that  we  choose  the  parameters  such  that  the  maximum 
fn  at  an  interface  is  of  the  same  magnitude  as  that  observed 
in  relevant  uniaxial  tensile  experiments  (Vratsanos  and  Farris, 
1993).  The  specific  parameter  values  chosen  are  E„  =  10^  MPa, 
E(  =  10"^  MPa  and  Au^  =  lO'^  m  unless  specified  otherwise, 

A  4-node  bilinear  hybrid  element  with  two  effective  degrees 
of  freedoms  (w,  v)  is  used,  along  with  the  interface  element 
described  in  Section  3.2.  The  numerical  integration  of  the  inter¬ 
face  element  stiffness  matrix  is  carried  out  by  a  2-point  Gaussian 
quadrature.  Due  to  the  instability  of  the  debonding  process,  all 
the  numerical  examples  are  displacement  controlled. 

To  validate  the  numerical  scheme  and  to  gain  insight  into 
the  behavior  of  composites  undergoing  the  evolution  of  micro¬ 
damage,  we  analyze  several  problems.  First,  computations  are 
earned  out  for  unit-cell  problems  to  evaluate  the  influence  of 
bond  quality,  the  effect  of  the  inclusion  size  (volume  fraction), 
shape  and  orientation,  as  well  as  the  boundary  loading  and 
constraints.  Second,  the  interaction  between  multiple  inclusions 
of  identical  or  different  sizes  is  investigated. 

In  the  following,  stress  refers  to  Cauchy  stress  if  not  specified 
otherwise,  strain  connotes  nominal  strain. 

4.1  Uniaxial  Stretching  of  an  Elastomer  Containing  a 
Circular  Inclusion.  A  3  mm  X  3  mm  square  elastomer  ele¬ 
ment  containing  a  circular  inclusion  (radius  1.3  mm,  volume 
fraction  59  percent)  is  analyzed  in  a  sequence  of  displacement 
controlled  deformations.  In  all  cases,  the  boundary  condition  is 
mixed,  i.e.,  a  constant  normal  displacement  is  prescribed  on 
two  opposite  sides  and  the  remaining  two  sides  are  traction  free. 

First  the  “global”  stress-strain  behavior  of  the  element  is 
compared  for  different  bond  qualities  at  the  interface,  namely 
unbonded,  perfectly  bonded  and  bonded  according  to  the  in¬ 
terfacial  relation  (3).  These  overall  stress-strain  responses  are 


also  compared  with  that  of  a  pure  elastomer  square  and  that  of 
an  elastomer  square  with  a  hole  of  radius  1.3  mm  at  the  center. 

Figure  2  illustrates  the  global  stress-strain  response  of  a  single 
element  with  and  without  an  inclusion.  As  expected  one  notes 
that  compared  to  pure  elastomer,  the  configuration  with  a  per¬ 
fectly  bonded  inclusion  is  stiffer  and  that  with  a  hole  is  softer 
than  the  pure  elastomer.  In  addition  one  observes  that  the  ele¬ 
ment  with  a  bonded  inclusion  loses  stiffness  gradually  as  the 
debond  increases,  and  becomes  more  compliant  than  the  pure 
elastomer.  Also  one  observes  that  when  a  large  area  is  de- 
bonded,  the  geometry  behaves  like  one  with  an  unbonded  inclu¬ 
sion,  The  largest  nominal  strain  is  in  the  order  of  600  percent. 

4.2  Effect  of  Inclusion  Size.  The  effect  of  inclusion  size 
in  a  unit-cell  is  considered  next  and  documented  in  Fig.  3  for 
inclusion  radii  of  r  =  0.25, 0.5,  1.0, 1.3  mm  in  an  initial  domain 
of  3  mm  X  3  mm  (Corresponding  volume  fractions  are  2.2, 
8.7,  35.0,  59.0  percent).  When  the  size  of  inclusion  is  small 
(0.25  mm,  0.5  mm,  1  mm)  relative  to  the  domain  size,  the 
overall  stress-strain  curve  is  monotonic;  when  the  size  of  the 
inclusion  is  large  (1.3  mm),  the  overall  stress-strain  trace  is 
not.  The  larger  the  inclusion  size,  the  sooner  the  specimen  de¬ 
bonds  and  loses  stiffness,  because  debond  initiation  occurs  ear¬ 
lier  for  the  larger  inclusion.  When  the  inclusion  size  is  very 
small  (0.25  nun),  the  stiffness  of  the  composite  element  is 
larger  than  that  for  pure  elastomer  even  when  debonding  occurs. 

4.3  Change  of  Boundary  Constraints.  To  check  the  ef¬ 
fect  of  boundary  conditions  on  the  overall  behavior  of  a  unit¬ 
cell  (volume  fraction  35  percent),  we  change  the  boundary 
conditions  at  the  “sides”  of  the  unit  cell  from  traction  free,  to 
(a)  elastic  foundation  (i.e.,  the  side  walls  of  the  unit  cell  are 
constrained  horizontally  by  springs),  (b)  zero  normal  displace¬ 
ment  (wi  =  0).  The  behavior  of  the  unit-cell  changes  noticeably 
(Fig.  4) :  over  the  strain  range  presented,  the  traction  free  condi¬ 
tion  leaves  part  of  the  bond  between  the  matrix  and  inclusion 
intact,  while  the  lateral  displacement  constraint  leads  to  early 
dewetting  closely  followed  by  the  elastic  constraint,  but  with 
somewhat  later  unbonding. 

4.4  Effect  of  Inclusion  Shape  and  Orientation.  Inclu¬ 
sions  are  generally  not  round.  The  effect  of  the  shape  of  an 
inclusion  is  investigated  by  considering  inclusions  in  the  shape 
of  a  circle,  an  ellipse  and  a  square  located  in  the  center  of  a 
square  elastomer  matrix.  The  shape  of  an  inclusion  affects  the 
time  of  debond  initiation  and  the  growth  of  voids  (detail  results 
not  reported  here). 

The  orientation  effect  is  investigated  by  considering  a  speci¬ 
men  with  an  elliptic  inclusion  (aspect  ratio  2,  volume  fraction 


Fig.  2  Evaluation  of  bond  qualities:  overall  stress-strain  curves 
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17.5  percent)  in  three  orientations  (nnajor  axis  is  oriented  90, 
45  and  0  degrees  with  respect  to  the  tension  axis).  Figure  5 
shows  that  when  the  elliptic  inclusion  is  oriented  at  90  or  45 
degrees,  the  overall  stress-strain  curves  are  non-monotone.  For 
the  90  deg  case,  the  debonding  occurs  over  a  large  portion  of 
the  interface  almost  at  once,  which  corresponds  to  the  drop  in 
stress  in  Fig.  5  (curve  3).  A  similar  situation  arises  for  the  45 
deg  case.  These  indicate  that  when  the  largest  cross  section  of 
the  inclusion  is  perpendicular  to  the  loading  axis,  the  debonding 
occurs  earlier  and  faster  than  for  (near)  alignment. 

4.5  Interaction  of  Inclusions.  A  comparison  was  also 
made  between  different  numbers  of  inclusions  in  a  block  of 
elastomer  which  have  the  same  volume  fraction  and  identical 
inclusions.  The  number  of  inclusions  considered  are  1,2,  and 
4.  It  is  observed  that  there  is  basically  no  difference  between 
the  case  of  one  inclusion  and  that  with  two  inclusions.  The 
overall  stress-strain  curve  for  the  elastomer  with  four  inclusions 
differs  little  from  those  with  one  or  two  inclusions,  and  there 
is  a  trend  that  the  overall  stress-strain  curve  converges  to  a 
single  curve  as  the  number  of  inclusions  increases. 

A  square  elastomer  with  4  and  16  inclusions  was  also  consid¬ 
ered  for  which  the  inclusion  sizes  fluctuate  around  a  mean  value 
(0.3  mm  for  the  4-inclusion  case,  0.4  mm  for  the  16-inclusion 
case).  Again,  normal  displacements  are  described  at  two  oppo¬ 
site  sides  of  the  elastomer.  The  remaining  two  sides  are  traction 


Fig.  4  Effect  of  changing  boundary  constraints 


(a) 

Fig.  5(a) 


(b) 

Fig.  5(b) 

Fig.  5  Effect  of  orientations:  (a)  the  overall  stress-strain  curves;  (b)  a 
stress  contour 


free.  The  size  of  the  elastomer  is  2  mm  X  2  mm  for  the  4- 
inclusion  case,  and  4  mm  X  4  mm  for  the  16-inclusion  case. 
For  the  4-inclusion  case,  the  computation  extends  to  18.3  per¬ 
cent  overall  strain  (the  largest  normal  strain  in  the  elastomer  is 
43  percent),  after  which  it  diverges.  For  the  16-inclusion  case 
the  computation  reaches  4  percent  overall  strain  with  the  largest 
normal  strain  being  23  percent,  thereafter  the  algorithm  diverges 
(see  Fig.  6). 

A  more  detailed  examination  of  the  results  for  the  two  cases 
shows  that  debonding  occurs  first  at  the  larger  inclusions  and 
voids  formed  there  grow  larger  as  the  external  load  increases 
compared  to  those  at  smaller  inclusions.  The  voids  formed  at 
smaller  inclusions  seem  not  to  grow  or  grow  very  slowly.  This 
is  consistent  with  the  observation  made  in  the  unit-cell  situation. 
What  cannot  be  observed  in  a  unit-cell  setting  is  that  the  voids 
formed  at  the  inclusion-matrix  interface  can  grow,  they  can  also 
shrink  as  external  loading  is  increased  due  to  the  interaction 
among  inclusions  (see  Fig.  7),  and  that  a  sufficiently  large  void 
tends  to  coalesce  with  a  sufficiently  large  nearby  void,  but  not 
necessarily  the  nearest  void.  In  addition  discrete  segmental  in¬ 
creases  are  observed  in  the  overall  stress-strain  curve  for  the 
16-inclusion  case,  but  not  for  the  4-inclusion  case,  this  is  due 
to  the  fact  that  the  16-inclusion  case  has  a  larger  volume  fraction 
of  inclusions  than  the  4-inclusion  case  (about  50  percent  versus 
30  percent).  This  “roughness”  in  the  stress-strain  response  has 
also  been  observed  by  Gonzalez  and  Knauss  ( 1997),  where  the 
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volume  fraction  of  (three-dimensional)  particles  is  about  70 
percent. 


5  Some  Numerical  Issues 

It  is  observed  that  the  convergence  of  computations  is  very 
sensitive  to  the  mesh,  time  increment  and  quadrature  for  the 
integration  of  the  interface  element  stiffness  matrix.  When  de¬ 
bonding  occurs,  a  very  large  number  of  iterations  is  necessary. 
We  tested  Gaussian  quadrature  and  Lobatto  quadrature  for  the 
integration  of  the  interface  element.  In  the  present,  large  defor¬ 
mation  occasion,  Gaussian  quadrature  was  found  more  satisfac¬ 
tory  than  Lobatto  quadrature,  which  is  more  efficient  in  cases 
of  small  matrix  deformation  (see,  for  example,  Schellekens  and 
De  Borst,  1993). 

In  order  to  analyze  a  force-controlled  problem  setting  or  an 
unsuble  failure  process  in  a  displacement  controlled  problem 
setting,  a  dynamic  rather  than  a  quasi-static  formulation  is 
needed.  The  interfacial  decohesion  propagates  very  rapidly  dur¬ 
ing  an  unstable  failure  process  as  observed  by  Gonzalez  and 
Knauss  (1997).  The  unstable  nature  of  a  fast  propagating  de¬ 
cohesion  also  shows  up  in  stress-strain  curves  as  a  small,  sudden 
drop  in  the  stress. 


Fig.  6(a) 


(b) 

Fig.  6{b) 

Fig.  6  Interaction  of  (16)  inclusions  of  different  sizes:  (a)  the  overall 
stress-strain  curve;  (h)  a  stress  contour 


Fig.  7  A  4-inclusion  case:  voids  formed  at  interfaces  can  grow  and/or 
shrink  when  external  load  increases  (a)  overall  strain  8.6  percent;  (b) 
overall  strain  18.3  percent 


6  Conclusions 

By  characterizing  the  particles,  matrix  and  interfaces  in  a 
particulate-filled  elastomer  individually,  the  whole  failure  pro¬ 
cess  (interfacial  debonding,  growth  of  voids  formed  at  inter¬ 
faces)  is  analyzed  by  the  FEM.  It  is  demonstrated  that  the 
scheme  can  predict  the  damage  evolution  in  particulate-filled 
elastomer  in  a  “discontinuum”  way.  The  predictions  are  quali¬ 
tatively  the  same  as  those  observed  in  laboratory  experiments. 

The  main  observations  are  that  the  behavior  of  particulate- 
filled  elastomers  depends  on  the  size,  shape,  orientation  and 
interface  properties:  in  an  elastomer  containing  regularly  distrib¬ 
uted  identical  particles,  the  overall  stress-strain  curve  seems  to 
converge  to  a  single  curve  as  the  number  of  particles  increases 
(fixed  volume  fraction);  in  a  filled-elastomer  containing  differ¬ 
ent  size  particles,  the  interaction  among  particles  plays  an  im¬ 
portant  role  in  the  growth  and  coalescence  of  voids.  The  interfa¬ 
cial  failure  process  is  highly  unstable,  especially  when  the  vol¬ 
ume  fraction  of  inclusions  is  relatively  large. 

In  order  to  analyze  more  complicated  problems,  such  as  inter¬ 
action  of  particles  and  a  macroscopic  crack,  propagation  of  a 
crack  via  the  growth  and  coalescence  of  voids,  the  computa¬ 
tional  approach  must  be  made  more  robust  and  stable. 
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SUMMARY 

The  stability  of  homogeneous  phase  separation  in  finite  solids  containing  planar, 
cylindrical  or  spherical  interfaces  is  investigated  analytically.  Explicit  stability  conditions 
are  deduced  for  each  interface  geometry.  It  is  shown  how  the  interaction  of  load  (force  or 
displacement),  material  properties  of  the  phases  and  interface  properties  jointly 
determine  the  stability  of  the  interface  separation  process. 

INTRODUCTION 

In  an  effort  to  model  failure  and  fracture  in  highly  particle  filled  solids  (solid  propellant 
rocket  fuels  for  example)  it  has  been  proposed  to  characterize  the  damage  in  the  crack  tip 
region  through  discretized  representation  of  the  deformations  around  particles,  including 
their  separation  from  binder,  which  is  usually  refereed  to  as  "dewetting"  (Zhong  and 
Knauss  [1]).  The  extent  of  the  region  in  which  this  phenomenon  needs  to  be  modeled  in 
detail  is  determined  experimentally  to  be  on  the  order  of  a  millimeter  or  a  fraction 
thereof.  Outside  this  domain  the  composite  is  to  be  described  by  a  (possibly  nonlinear) 
continuum  constitutive  law.  Thus  the  individual  phases  of  the  composites  are 
characterized  separately,  i.e.  matrix  material,  particles  and  the  interfaces  are  assigned 
different  constitutive  behavior  in  a  detailed  finite  element  model. 

In  evaluating  such  models  for  various  phase  properties  and  volume  fractions,  instabilities 
occur  that  raise  the  question  whether  their  origins  are  numerical  or  physical  in  nature. 
Since  interfaces  are  ubiquitous  in  (fibrous,  particulate  or  layered)  composites  and  damage 
accrued  through  interfacial  failure  is  an  important  phenomenon  in  these  solids,  we 
address  the  issue  how  the  interaction  of  specimen  size,  filler  geometry,  properties  of 
matrix  and  interfacial  constitution  lead  to  stable  vis-a-vis  imstable  interfacial  separation. 
It  thus  the  purpose  of  this  note  to  illuminate  this  question  through  deriving  global  criteria 
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that  establish  conservative  estimates  for  the  onset  of  instability.  In  the  context  of  layered, 
fibrous  and  particulate  composites  we  consider  thus  finite  domains  containing  planar, 
cylindrical  or  spherical  interfaces. 

Interfaces  can  be  characterized  at  least  qualitatively  by  interfacial  constitutive  models  ( 
see  for  example,  Needleman  [2],  Ungsuwarrungsri  and  Knauss  [3]).  Suo,  Ortize  and 
Needleman  [4]  investigated  the  stability  of  solids  containing  interfaces  in  terms  of  two 
semi-infinite  solids  bonded  along  a  plane  by  defining  the  onset  of  instability  as  the 
appearance  of  non-unique  solutions  to  the  associated  boundary-value-problem.  In  the 
sequel,  the  principle  of  minimum  potential  energy  is  used  to  select  a  favorable  solution 
when  bifurcation  occurs.  The  unstable  solution  is  considered  in  the  Liapunov  sense  such 
that  an  infinitesimal  change  in  the  excitation  produces  a  finite  change  in  response.  It  is 
shown  that  the  onset  of  non-uniqueness  of  solutions  to  boundary  values  problems  does 
not  necessarily  lead  to  instability  in  an  energy  sense,  a  result  that  includes  the  special 
situation  of  a  propagation  craze-crack  employing  more  specific  interfacial  constitutive 
descriptions  (Ungsuwarungsri  and  Knauss  [5]). 

Stability  analyses  for  the  development  of  inhomogeneous  deformation  fields  have  been 
performed  in  connection  with  shear  band  formation  when  coupled  with  material  strain 
softening  behavior  as  (in  the  absence  of  any  fracture),  for  example,  by  Needleman  [6]. 
More  recently  Levy  [7]  has  studied  how  phase  concentration  affects  stability  through  a 
bifurcation  condition  in  a  two  dimensional  setting  and  for  a  specific  interface  model.  Here 
we  consider  a  more  general  model  of  interface  cohesion,  and  identify  how  the  interaction 
of  load  (displacement  or  force),  matrix  material  properties  and  interface  properties  jointly 
determine  the  stability  of  an  interfacial  separation  process.  We  show  specifically  that 
interfacial  separation  instability  involves  the  maximum  negative  slope  of  the  unloading 
portion  of  the  interface  cohesion  model  together  with  the  matrix  modulus  and  the 
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specimen  size.  Thus  the  work  reported  here  provides  new  perspectives  for  numerical 
simulation  of  interfacial  failure  process. 


AN  INTERFACIAL  CONSTITUTIVE  MODEL 

Interfacial  constitutive  models  considered  to  date  are  typically  rate  independent  and  relate 
the  separation  to  the  traction  across  an  interface.  The  general  features  of  these  models  are 
that  the  traction  across  the  interface  first  increases  and  then  decreases  with  interfacial 
separation;  by  definition  total  separation  at  the  interface  is  achieved  when  the  traction 
across  the  interface  has  dropped  to  zero.  For  bond-normal  separation,  an  interfacial 

constitutive  relation  between  the  traction  T  and  separation  A 

r  =  /(A)  (1) 

is  illustrated  schematically  in  Figure  1. 

A  PLANAR  INTERFACE 

Consider  an  elastic  two-phase  solid  consisting  of  two  rectangular  plates  (with  Young’s 
modulus  and  Poisson’s  ratios  E, ,  v, ,  Ej ,  Vj )  bonded  as  sho'ft'n  in  Figure  2.  The  interface 
is  vanishing  thin  and  is  characterized  by  the  interfacial  constitutive  relation  of  the  form 
described  in  Figure  1.  Because  the  system  obviously  becomes  unstable  as  soon  as 

T  =  ,  if  normal  traction  t_  are  prescribed  at  X2  =  ±L.  We  consider  here  only  a 

displacement  control  problem  with  boundary  conditions 
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/  =  0, 


at  X,  =  0,a , 


=0,  M2  =+M,  at  X2=±L. 

For  reasons  of  simplicity  eliminate  shear  stress  by  letting  v,  =  Oand  ^3=  0.  Equilibrium 

(plane  stress)  and  compatibility  determine  that  the  edge  displacement  u,  the 
homogeneous  stress  <j{=(72t)  and  the  separation  A  at  the  interface  are  related  by 

/(A)  =  ^(2m-A)  (2) 

where  E  =  / (£,  +  £2) • 

This  nonlinear  algebraic  equation  can  be  analyzed  geometrically.  Consider  the 

intersection  of  the  curve  y  =  /(A)  and  the  straight  line(s)  :v(m,  A)  =  2E(2u-A)IL.  The 
number  of  solutions  to  (2)  (intersections  of  the  two  y-functions)  depends  on  the  slope  of 

the  straight  line(s)  and  the  magnitude  of  u  for  given  interface  properties  (See  Figure  3)  : 

When  the  slope  {2E  /  L)  of  straight  line(s)  is  sufficiently  large,  the  solution  to  (2)  is 

unique  for  all  u ;  for  smaller  slopes,  the  solution  can  bifurcate  when  u  is  large.  Thus  the 
condition  for  the  existence  of  a  unique  solution  to  (2)  is  thus 

9  F 

—  >  max(-/-(A)).  (3) 

0<A<Ao 

Let  n  be  the  total  energy  (work  done  on)  in  the  system,  consisting  of  the  sum  of  the 
energy  stored  elastically  in  the  plates  and  the  work  done  against  the  interface  traction 
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(4) 


n = + ]mdi, 

a  2E  I  2.E  0 

with  /  given  by  (2). 


^  >  jnax  (-/'(A)) ,  the  total  energy  n  of  the  system  is  a  monotonic  function  of 

L  0<A<Ao 

the  applied  displacement  u  ;  When  ^  sufficiently  large, 

there  are  two  or  three  solutions  to  (2).  For  each  solution,  the  corresponding  n  is  a 
monotonic  function.  The  acceptable  solution  is  that  which  minimize  the  potential  energy. 


It  is  thus  seen  that  when  £  / 1  is  large,  a  stable  homogenous  decohesion  of  the  interface 
is  possible  in  the  sense  that  a  continuously  varying  separation  of  the  interface  occurs  as  a 
function  of  the  applied  displacement;  otherwise  the  instability  of  the  system  occurs.  This 

can  happen  even  before  the  stress  at  the  interface  reaches  and  the  two  plates  separate 

suddenly  at  a  non  zero  stress,  provided  1  / 1  is  sufficiently  small.  In  other  words,  the 
stiffer  the  bulk  materials  and  smaller  the  domain  size,  the  more  stable  the  system  is  under 
boundary  displacement  control. 


A  CIRCULAR  INTERFACE 

Consider  next  a  simple  case  of  a  cylindrical  rod  or  fiber  (assumed  rigid)  embedded  in  an 
elastic  cylindrical  matrix  with  inner  radius  ci ,  outer  radius  b ,  possessing  an  elastic 
modulus  E  and  Poisson’s  ratio  v  under  plane  strain  conditions,  assuming  that  the  rod 
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(inclusion)  is  rigid.  The  interface  is  again  characterized  by  the  interfacial  constitutive 
relation  illustrated  in  Figure  1. 

Assuming  axisymmetric  deformation  there  are  no  shear  stresses  in  the  system. 

The  boundary  conditions  (polar  coordinates  (r,6>))  are  then 

(I)  r  =  b,u^=u,t^  =  0.  (displacement control) 
or 

(II)  r  =  b,t^=~t,t0=O.  (force  control) 


It  is  easy  to  show  that  the  stability  conditions  for  these  two  boundary  conditions  are 
(1)  Displacement  control 

2Eb(b^-a^)  b\\- v)  +  a\l+ v)  ^  ^ 

- - — - - >  max(-  /  (A) 

(a^ +6^)^(1- v^)  o<a<a„ 


(5) 


(2)Traction  control 


^->  n.ax(-/'( 

aCaHl-vj+i-Hl  +  v))  ’  0<A<Ao^ 


(6) 


The  potential  energy  as  a  function  of  applied  displacement  is  given  as  an  example  in 
Figure  4  for  an  unstable  decohesion  process.  Compared  with  the  case  for  the  planar 
interface,  we  see  that  (1)  stable  interfacial  decohesion  is  possible  even  the  process  is  force 
controlled,  (2)  debonding  process  is  more  stable  the  larger  b/a  is;  this  is  opposite  to  the 


7 


finding  for  the  planar  interface  where  a  smaller  “L”  corresponds  to  the  more  stable 
debonding  process.  Thus  the  stability  of  interfacial  failure  (dewetting)  depends  on  the 
interfacial  properties  as  well  as  on  the  geometry  of  the  interface  and  loading  conditions. 


A  SPHERICAL  INTERFACE 

A  three  dimensional  generalization  of  the  results  for  the  circular  interface  case  is  easily 
realized  by  considering  a  rigid  spherical  particle  (r  =  a)  embedded  in  a  spherical  domain 
(bonded  by  r  =  b  )  of  Young’s  modulus  E  and  Poisson’s  ratio  v  subjected  to 
axisymmetric  deformation.  This  consideration  involves  a  spherical  interface,  the 
geometrically  simplest  interface  in  particulate  composites. 

The  boundary  conditions  (spherical  coordinates  (r,6,<p))  are  then 


(I)  r  =  b,u^  =u,tg  =  0,t^  =  0 .  (displacement  control) 


or 


(II)  r  =  b,t,=  t,tg  =  =  0 .  (force  control) 


leading  to  the  stability  conditions 


(I)  Displacement  control 


2(1  2v)6^+(l+v)a  E{b — O  >  j^ax  (-/’(A)) . 

a\2(2-v)b^ -i}  +  v)a^)  b-a  0<A<Ao 


(7) 
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(II)  Traction  control 


E{b^  -a^) 

a{{\  +  v)b^  +(\-2v)a^) 


>  max  (-/'(A)). 

0<A<Ao 


(8) 


DISCUSSION 

The  above  stability  conditions  are  consistent  with  recent  numerical  experiences  (Zhong 
and  Knauss  [1],  Needleman  [8])  for  both  planar  and  circular  interfaces  in  much  more 
complex  situations.  As  far  as  the  stability  of  a  system  involving  homogeneous  interfacial 
decohesion  is  concerned,  there  is  not  much  difference  between  a  circular  interface  and  a 
spherical  interface. 

The  instability  is  associated  to  the  sudden  drop  of  traction  at  interfaces  and  associated 
energy  release  when  interfacial  failure  occurs.  This  local  instability  poses  considerable 
difficulties  for  the  numerical  simulation  of  interfacial  failure  in  a  quasi-static  setting  and 
will  lead  to  the  divergence  of  numerical  computations.  One  way  to  cope  with  this 
complication  is  to  model  the  interfacial  failure  process  as  a  dynamic  process.  Another 
way  is  to  “regularize”  the  interfacial  constitutive  relation  at  the  point  of  the  instability 
onset,  replacing  the  original  traction  separation  curve  by  a  “modifying”  line  derived  from 
stability  analysis.  Such  a  regularization  may  lead  to  avoidance  of  instabilities  in 
numerical  computation  while  retaining  the  basic  feature  of  the  original  response.  This 
regularization  makes  the  corresponding  numerical  simulation  somewhat  similar  to  the 
“nodal  release”  technique'. 


*  In  this  regularization  process  the  fracture  energy  remains  as  the  same  as  that  defined  by  the  interfacial 
model  while  the  fracture  energy  associated  with  the  “nodal  release”  technique  is  arbitrary  (undefined). 
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We  have  considered  here  only  uniform  decohesion  at  interfaces.  A  more  complete 
description  would  concern  the  stability  of  the  local  failure  process  through  an 
inhomogeneous  decohesion.  Intuitively,  an  interfacial  failure  process  involving 
inhomogeneous  decohesion  is  more  stable  than  that  involving  homogeneous  decohesion 
here,  because  the  surrounding  bonded  region  tends  to  stabilize  the  decohesion  process. 
Thus  the  stability  estimation  obtained  here  for  homogenous  decohesion  might  provide  a 
non-conservative  estimate  of  the  onset  of  instability  for  an  arbitrary  interfacial  failure 
process.  It  has  the  advantage  of  relating  in  a  simple  manner  the  domain  size  and  the 
properties  to  a  stability  criterion,  a  situation  that  requires  much  more  analysis  for  an 
inhomogeneous  dewetting  process. 
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EFFECTS  OF  PARTICLE  INTERACTION  AND  SIZE  VARIATION 
ON  DAMAGE  EVOLUTION  IN  FILLED  ELASTOMERS 
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Pasadena,  CA  91125,  USA 

Abstract 

A  micromechanical  analysis  of  damage  evolution  (interfacial  debonding)  in  particle- 
filled  elastomers  addresses  the  effect  of  the  interactions  between  particles  and  of  vari¬ 
ation  in  filler  size.  The  composite  is  treated  as  an  assembly  of  two  constituents  in  a 
finite  element  model.  It  is  shown  that  the  interaction  between  particles  controls  the 
damage  evolution:  (1)  For  high  volume  fraction,  a  relatively  small  change  in  particle 
size  has  a  surprisingly  large  effect  on  local  material  response;  (2)  for  large  differences 
in  particle  sizes(e.g.  bimodal  distribution),  damage  occurs  at  interfaces  between  large 
particles  and  matrix,  with  limited  damage  occuring  at  small  particles.  While  these 
effects  of  particle  interaction  and  size  variation  tire  smoothed  out  in  a  large  ensemble 
of  particles,  it  is  foreseeable  that  they  are  an  important  factor  in  a  failure  process  such 
as  maccroscopic  crack  propagation,  which  spans  scales  considerably  larger  than  the 
mtiximum  particle  size.  Specifically,  one  expects  thus  that  in  the  vincity  of  a  macro¬ 
scopic  crack  the  large  particles  become  sites  for  small  cracks  which  coalese  to  larger 
ones  and  join  up  with  the  macro  crack,  while  small  particles  operate  primarily  so  as  to 
locally  stiffen  the  matrix  without  incurring  significant  damages  in  their  vincity. 


1  Introduction 

Recently  substantial  efforts  have  been  made  in  community  of  mechanics  of  materials  towards 
damage  evolution  in  composites  and  its  effect  on  global  material  behavior.  Theoretical  stud- 
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ies  with  constitutive  models  accounting  for  damage  have  been  proposed;  see  for  example, 
Farris  and  Schapery  (1973),  Schapery(1986),  Govindjee  and  Simo  (1992),  Vratsanos  and 
Farris  (1993),  and  Ravichandran  and  Liu  (1995).  By  proposing  a  constitutive  model,  one 
implicitly  treats  a  damaged  material  as  a  homogenized  material  and  neglects  any  detailed 
interaction  between  various  damage  sites.  In  attempts  to  inquire  into  the  mechanics  of  the 
damage  development  finite  element  analyses  tend  to  make  use  of  unit-cell  models  (Needle- 
man  (1987),  Yeh  (1992),  and  Walter,  Ravichandran  and  Ortiz  (1997)  in  which  one  assumes 
that  fillers  (fibers,  particles  etc)  are  uniformly  of  the  same  size  and  (periodically)  distributed 
homogeneously  in  the  composite.  However  the  unit-cell  approach  cannot  address  the  inter¬ 
action  between  particles  or  their  clustering  of  fillers,  nor  can  it  be  instructive  on  the  effect 
of  variations  in  filler  size. 

The  motivation  behind  the  current  work  is  the  desire  to  develop  a  computational  scheme 
to  deal  with  the  fracture  behavior  of  particulate  composites.  To  this  end  Zhong  and  Knauss 
(1997b)  have  proposed  a  hybrid  discrete-continuum  approach  which  incorporates  the  discrete 
damage  evolution  in  filled  elastomers.  The  basic  thought  is  to  characterize  the  material  far 
from  the  crack  tip  as  a  homogeneous  body,  while  close  to  the  crack  tip,  in  a  region  measured 
possibly  in  millimeters  or  thereof  the  discrete  particulate  interaction  with  damage  evolution 
are  modeled  as  illustrated  in  Figure  1.  Such  pursuits  appear  feasible  in  light  of  the  ever  more 
rapidly  increasing  power  of  computing  machines.  Initial  studies  directed  to  that  goal  have 
shown  that  for  small  numbers  of  particles  the  material  behavior  is  reproduced  qualitatively 
(Zhong  and  Knauss  (1997b)). 

In  this  paper  we  concentrate  on  the  effect  which  the  interactions  of  particles  have  on  the 
damage  evolution  along  with  the  quest  for  understanding  the  role  which  particle  size  plays 
in  this  picture.  For  the  purpose  of  this  initial  study,  the  elastomer  matrix  is  characterized  as 
a  hyperelastic  Ogden  material  (Ogden,  1972),  although,  in  reality,  the  elastomer  may  posess 
time  dependent  characteristics  Because  the  elastomer  matrix  is  typically  soft  relative  to 
the  particles,  the  latter  may  be  considered  to  be  rigid  The  interface  between  particle  and 

^The  motivation  for  this  work  comes  from  the  failure  behavior  of  solid  propellant  rocket  fuel;  these 
materials  possess  time  or  rate  dependent  properties,  but  for  exploratory  purposes,  it  suffices  to  ignore  the 
effects. 

^This  choice  is  made  for  convenience.  Under  high  rate  loading,  the  binder  may  be  stiff  by  compression 
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matrix  can  fail  and  is,  therefor,  modeled  by  a  cohesive-zone  type  feature  with  the  bond  at  the 
interfaces  so  ’’weak”  that  only  interfacial  debonding  needs  to  be  considered  (  no  cavitation 
occurs  in  the  elastomer).  A  two  dimensional  version  of  the  problem  is  considered  here  and  a 
fully  three  dimensional  problem  description  is  left  to  future  investigations,  since  it  requires 
only  a  larger/faster  computer. 

In  the  sequel  we  delineate  the  interface  model  along  with  the  associated  stability  issue  in 
Section  2.  The  computational  model  is  presented  in  Section  3  including  the  FEM  formulation 
for  hyperelasticity  and  the  interface  element.  To  evaluate  the  effect  of  the  interaction  between 
particles  (inclusions)  of  different  size,  we  analyze  configurations  of  four  (4)  as  well  as  sixteen 
(16)  inclusions  wherein  the  particle  sizes  vary  by  small  amounts  (10%).  In  addition  we 
examine  the  interaction  between  inclusions  for  a  bimodal  size  distribution  and  find  that  (1) 
for  a  high  volume  fraction  of  filler,  a  small  change  in  particle  size  has  a  large  effect  on  material 
response,  (2)  when  there  is  a  big  difference  in  particle  sizes(such  as  bimodal  distribution), 
damage  occurs  along  interfaces  of  large  particles,  little  damage  taking  place  around  small 
particles. 

2  An  interfacial  constitutive  model 

To  model  the  interfacial  separation,  we  imitate  the  cohesive  failure  by  a  traction-displacement 
relation  that  replaces  the  local  consititutive  behavior  in  a  thin  layer  (zero  thickness  in  the 
model).  The  interfacial  traction  is  assumed  to  increase  with  increasing  interfacial  separation 
to  a  critical  value,  after  which  the  traction  decreases  with  increasing  separation,  eventually 
vanishing,  which  condition  then  corresponds  to  total  separation,  (cf  Ungsuwarungsri  and 
Knauss;1987). 

2.1  The  interfacial  constitutive  model 

Let  n  be  the  normal  to  an  interface,  and  s,  t  the  corresponding  orthogonal  tangential  direc¬ 
tions;  En,  Eg  and  Et  are  elastic  moduli  of  a  nonlinear  (dissipative)  spring  in  the  directions 

SO  that  the  particles  will  also  deform.  That  case  is  computationally  treatable  as  well,  but  requires  larger 
computational  resources. 
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(n,s,t).  The  spring  response  is 


0  0 

D=  0  0  ,  (1) 

0  0 

as  long  as  the  normal  separation  Au„  is  less  than  the  critical  normal  separation  Auc-  Once 
Auji  >  Auc,  the  spring  moduli  degrade  such  that 

D=  0  0  (2) 

_0  0 

for  Auc  <  Aun  <  2Auc,  and  vanish  for  Au„  >  2Auc. 

The  nonlinear  interface  constitutive  relation  is  thus 

f  =  DAu(„,),  (3) 

with  Au  the  relative  displacement  in  local  coordinates  (n,  s,  t)  and  f  the  conjugate  interfacial 
traction  defined  on  the  undeformed  reference  configuration. 

Only  failure  due  to  normal  separations  is  accounted  for  in  the  current  model.  While 
one  needs  to  consider  shear  failure  at  an  interface  for  more  general  situations,  it  is  not 
difficult  to  include  shear  failure  into  the  above  interfacial  model  as  demostrated,  for  example 
by  Tvergaard  and  Hutchinson  (1993).  However,  for  the  type  of  solid  propellant  related 
problems  under  consideration,  the  tensile  strength  of  the  interface  is  much  lower  than  the 
shear  strength  so  that  considering  the  ’’opening  mode”  of  failure  alone  encompasses  the 
major  damage  contribution. 

2.2  Stability  of  an  interfacial  failure  process 

In  1992  Suo,  Ortiz  and  Needleman  presented  a  stability  analysis  for  solids  possessing  inter¬ 
faces  by  considering  two  semi-infinite  solids  bonded  along  a  plane.  Instability  of  the  system 
was  defined  as  the  onset  or  development  of  nonunique  solutions  to  the  associated  boundary- 
value-problem.  In  a  more  global  approach  Zhong  and  Knauss  (1997a)  recently  studied  the 
stability  of  finite  solids,  with  planar  or  circular  interfaces  under  homogeneous  decohesion 
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(i.e.  the  whole  interface  debonds  uniformly  at  the  same  time)  deriving  compact  stability 
criteria  for  a  general  form  of  the  traction-displacement  interfacial  relation 

/  =  /(A»).  (4) 

For  homogeneous  decohesion  at  a  bi-material  planar  interface  the  stability  condition  is 

(6) 

where  Aug  corresponds  to  total  separation  at  the  interface.  The  right  hand  side  of  inequality 
(5)  represents  the  maximum  value  of  when  Aw  is  in  (0,  Aug).  It  is  clear  that  the  size 
of  the  specimen  (L),  effective  elastic  modulus  {Ee/f  =  ExE^/^Ei  +  E2)),  where  Ei  and  E2 
are  the  elastic  moduli)  and  the  unloading  part  of  the  interfacial  relation(4)  determine  the 
stability  of  the  system.  Similar  stability  conditions  are  obtained  for  circular  and  spherical 
interface  (see  Zhong  and  Knauss  (1997a)).  For  the  interfaces  in  filled  elastomers  modeled 
here,  Ei  is  the  modulus  of  matrix  material,  E2  is  the  modulus  of  particles,  taken  here  to 
be  infinity  because  of  their  relatively  high  stiffness  compared  to  rubbery  bindder,  so  that 
Egff  =  El  (matrix  modulus).  It  is  thus  clear  from  the  above  stability  condition  (5),  that  the 
more  compliant  the  matrix  material  is,  the  less  stable  is  the  corresponding  composite  with 
respect  to  interfacial  failure.  Thus  particulate-filled  elastomers  tend  to  be  less  stable  than 
particulate  metallic  composites  with  respect  to  the  interfacial  failure  process. 

An  interfacial  failure  process  associated  with  inhomogeneous  decohesion  at  interfaces  is 
intuitively  more  stable  than  that  associated  with  a  homogeneous  decohesion.  This  is  so 
because  the  surrounding,  bonded  regions  tend  to  stabilize  the  decohesion  process.  Thus 
the  stability  conditions  obtained  by  Zhong  and  Knauss  (1997a),  though  applied  only  to 
homogeneous  decohesion  at  interfaces,  provide  non-conservative  estimates  for  the  onset  of 
instability  in  interfacial  failure  processes. 

The  instability  due  to  interfacial  decohesion  is  local  because  of  the  rapid  (sudden)  drop 
of  tractions  at  the  interfaces.  This  local  instability  poses  a  problem  for  the  numerical  simu¬ 
lation  in  a  quasi-static  setting  because  it  can  (and  does)  lead  to  the  divergence  of  numerical 
computations. 
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3  Formulation  of  the  computational  model 


We  allow  for  finite  deformations  and  use  a  three-term  Ogden  strain  energy  function  for  the 
constitutive  description  of  the  matrix  as  a  hyperelastic  material.  For  reference  purpose,  we 
summarize  here  the  hyperelasticity  theory  and  its  hybrid  FEM  formulation  briefly. 

3.1  Hyperelasticity  and  its  FEM  formulation 

Let  flo  be  the  undeformed  reference  configuration  of  the  hyperelastic  material (s),  Clt  the 
current  deformed  configuration,  and  Vo  and  Vt  the  corresponding  volumes.  XeVlo,  yeflt- 
If  u(X)  is  the  displacement  vector  at  material  point  X,  F  the  corresponding  deformation 
gradient  and  J  the  Jacobian  of  F  (unit  volume  change)  one  has 


=  X  +  u(X), 

(6) 

F=.^ 

ax’ 

(7) 

J  =  det{F). 

(8) 

Because  the  hyperelastic  material  is  (almost)  incompressible,  the  deformation  gradient 
F  is  multiplicatively  decomposed  into  its  deviatoric  and  dilatational  components,  F*^®"  and 
F’'®^  such  that 


p  _  pdevpvo/ 

•pdev  _ 

■pvol  ^  jl/3j^ 


(9) 

(10) 

(11) 


with  I  the  identity  tensor.  It  is  easy  to  check  that  <iet(F'^®’')  =  1,  det{F'’°‘)  =  J  so  that  F'^®" 
is  the  volume  preserving  part  of  F,  and  F"®*  its  dilatation  part. 

The  strain  energy  function  for  a  nearly  incompressible,  hyperelastic  solid  was  given  by 
(Ogden,  1984), 

W(F)  =  W(F‘^®’')  -h  (12) 


where  W  takes  account  of  the  deviatoric  deformation  and  (j)  with  (/>(!)  =  0,  of  dilational 
deformations. 
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Considering  frame  indifference,  and  further  assuming  that  the  hyperelastic  material  is 
isotropic,  one  has 

W{F)  =  uj{X,XAz)  +  <f>{J),  (13) 

where  Aj  are  the  principal  deviatoric  stretches. 

Ogden(1972)  assumes  that  u  has  the  form, 


u;(Ai,  A2,  A3)  =  E  ^(Ar  +  A?‘  +  A?‘  -  3)  (14) 

i=l 

so  that  the  small  strain  shear  modulus  is  /xq  =  IJ-i  and  the  relation  between  the  Cauchy 
stress  and  the  principal  stretch  is 


d\  '3^  '$X(  dJ  ’ 


(15) 


If  the  material  is  incompressible  (J  =  1),  one  has  (/>(J)  =  0,  and  the  corresponding 
stress-stretch  relation  is 

PiW 

(16) 


r.  -  X  -  P- 

1 1  ^  I 


dXi 


where  Pi  is  an  arbitrary  pressure  and  Aj  are  the  principal  stretches. 

For  nearly  incompressible  material,  the  usual  displacement  finite  element  formulation 
can  ’’behave  poorly”  because  the  effective  bulk  modulus  is  very  large  compared  to  its  shear 
modulus,  which  cause  the  stiffness  matrix  to  be  almost  singular  from  a  numerical  point  of 
view.  Consequently,  the  stress  calculated  at  the  numerical  integration  points  show  large 
oscillation  in  the  pressure.  To  avoid  this,  a  hybrid  element  or  mixed  formulation  has  been 
proposed  (see  for  example,  Simo  and  Taylor(1991)). 


3.2  Interface  element 

The  interfacial  constitutive  relation  is  implemented  by  constructing  an  interface  element  by 
the  principle  of  virtual  work  (Beer  (1985))  by  the  follwoing  steps;  (1) Set-up  local  coordinates 
(n,  s,  t)  to  describe  the  interface  geometrically;  (2)  Find  the  relative  displacements  at 
the  interface  in  local  coordinates;  (3)  Use  the  principle  of  virtual  work  and  the  interfacial 
constitutive  relation  to  obtain  an  equilibrium  condition  for  the  interface  element  in  terms  of 
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global  nodal  coordinates.  One  then  assembles  the  interface  element  stiffness  matrix  with  the 
remaining  element  stiffness  matrices  to  form  the  global  stiffness  matrix. 

If  F®  is  the  force  vector  for  nodes  of  the  interface  element,  a  is  the  element  displacement 
vector  and  B  the  transformation  matrix  between  nodal  displacements  and  displacement  in 
the  local  coordinates.  Then,  with  point  (3)  in  mind  one  has 


f  B^f  ds  =  F^ 

J  inter  face 


(17) 


or 


f  B'^DBads  =  F®, 

J  inter  face 


(18) 

f  interface 

If  n  refers  to  time  steps,  k  refers  to  the  sequence  of  iterations,  the  residual  force  due  to  the 
interface  element  is 

R  =  -  [  B^Dl^^Bdsai  (19) 

J  interface 

and  the  tangent  stiffness  matrix  of  the  interface  element  becomes 


K,  =  -f. 
*  du 


(20) 


For  the  interfacial  relation  (3),  the  tangent  stiffness  matrix  for  interface  elements  is  non- 
symmetric.  Using  the  tangent  stiffness  matrix  may  be  a  source  of  numerical  difficulty  since 
the  stiffness  matrix  for  the  bulk  material  is  symmetric.  We  use  direct  iteration,  which  is 
expressed  for  the  interface  elements  by 


/  B^D*Bds{Aa}  =  (21) 

J  inter  face  J  inter  face 

where  +  Aa.  The  implementation  of  the  FEM  formulation  for  hyperelasticity 

is  available  in  ABAQUS  5.5,  and  the  interface  element  (interfacial  constitutive  model;  the 
element  has  intially  zero  thickness)  is  incorporated  into  that  code  by  way  of  a  user  subroutine. 


4  Numerical  results 

We  consider  plane  strain  deformation  in  the  modeled  two-dimensional  composites.  The 
matrix  characteristics  are  determined  by  fitting  the  parameters  in  the  3-term  Ogden  function 
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to  the  experimental  data  of  Treloar(1940).  The  parameters  are  =  58.23MPa,  ai  =  1.029, 
fj,2  ==  1-35  X  10“^MPa,  0:2  =  10.707,  ^3  =  0.246MPa,  0:3  =  —2.957;  these  satisfy  the 
material  stability  condition.  The  parameters  in  the  interfacial  element  are  chosen  such  that 
the  maximum  fn  at  an  interface  is  of  the  same  magnitude  as  that  observed  in  relevant 
uniaxial  tensile  experiments  (Vratsanos  and  Farris(1993)).  The  specific  parameter  values 
chosen  are  En  =  10^ MPa,  Et  =  10‘^MPa  and  Auc  =  10“^m. 

The  numerical  integration  of  the  interface  element  stiffness  matrix  is  carried  out  by 
a  2-point  Gaussian  quadrature.  Due  to  the  inherent  proclivity  towards  instability  of  the 
debonding  process,  all  the  numerical  examples  are  examined  under  displacement  controll. 
As  a  convention,  we  mean  the  Cauchy  stress  when  writing  ’’stress”  unless  specified  differently 
and  the  word  ’’strain”  connotes  nominal  strain. 

The  common  approach  to  analyze  properties  of  composites,  whether  by  numerical  or 
closed  form  methods  involves  the  ’’unit  cell”  approach.  It  is  implicit  in  this  approach  that 
(a)  filler  particle  sizes  are  ’’distributed”  uniformly  and  periodically  and  (b)  the  unit  cell 
characteristics  are  adequatly  indicative  of  global  response  of  the  composite.  With  the  intent 
of  questioning  the  latter  assumption  we  examine  thus  particle  size  variations  and  divide  the 
effort  arbitrarily  into  consideration  of  small  and  large  variations  in  particle  size. 

4.1  Interaction  of  inclusions;  small  variations  in  inclusion  sizes 

Three  cases  are  considered:  composites  with  low  volume  fraction  (8.7%),  intermediate  volume 
fraction  (35.0%)  and  high  volume  fraction  (68.4%)  of  particles^.  In  these  three  cases  we 
analyze  a  square  domain  containing  4-inclusions  with  different  radii.  For  reference  purpose, 
we  also  analyzed  the  counterparts  of  these  problems  with  identical  inclusions,  the  radii  of 
which  are  the  mean  of  the  inclusion  radii  for  the  non-uniform  configuration.  In  the  seques  we 
refer  to  the  geometry  with  identical  particle  size  as  case  1  and  to  the  situation  with  different 
particle  size  as  case  2.  The  radii  of  inclusions  and  their  mean  values  for  each  case  are  listed 
in  table  1;  where  we  associate  the  case  number  1,  2,  3  with  volume  fractions;  particle  size 
variations  are  sub-cases  as  in  1.1  and  1.2  for  example  for  case  1. 

®The  upper  limit  of  volume  fraction  is  78.5%,  for  uniformly  distributed  ’’rods”  in  a  (two-dimensionl) 
composite  and  52.4%  for  uniformly  distributed  spherical  particles 
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Table  1:  Inclusion  radii  and  mean  radii 


Case 

Inclusionl 

Inclusion2 

Inclusions 

Inclusion4 

Mean 

volume 

radius 

radius 

radius 

radius 

radius 

fraction 

1 

0.49 

0.51 

0.52 

0.48 

0.5 

8.7% 

2 

1.1 

1.0 

0.99 

1.05 

1.035 

35.0% 

3 

1.41 

1.42 

1.39 

1.38 

1.4 

68.4% 

The  boundary  conditions  are  imposed  according  to  Figure  2:  Two  opposing  sides  of  a 
square  elastomer  are  displaced  uniformly  in  opposite  directions  by  a  prescribed  amount  and 
the  remaining,  latera  surfaces  remain  traction  free,  for  the  case  of  uniform  inclusions,  we 
take  advantages  of  symmetry  conditions. 

Low  volume  fraction  of  inclusions  (case  1) 

Comparing  stress  distribution  for  case  1.1  (inclusions  with  same  radius)  and  case  1.2 
(inclusions  with  differnt  radii)  at  several  global  strain  stages  (roughly  5%,  10%  and  16%), 
there  is  no  significant  difference  between  the  two  subcases.  This  observation  is  further 
supported  by  the  little  difference  between  global  stress  strain  response  of  these  two  subcases. 
(As  the  volume  fraction  is  low,  the  difference  between  the  particulate  composites  and  the 
pure  elastomer  is  small,  too).  These  results  were  reported  earlier  by  Zhong  and  Knauss 
(1997b). 

These  observations  indicate  that  when  the  volume  fraction  of  inclusions  of  a  composite 
is  small,  the  interaction  between  particles  is  weak,  and  thus  unit-cell  approach  or  damage- 
incorporating  constitutive  model  can  be  applied  with  confidence  even  though  the  sizes  of 
fillers  vary  a  little. 

Intermediate  volume  fraction  of  inclusions  (case  2) 

Case  2.1:  Inclusions  with  same  radius 

We  observe  that  1)  the  vacuoles  formed  at  the  interfaces  between  the  matrix  and  the 
four  inclusions  are  the  same  (obvious  due  to  the  symmetry  condition);  2)  the  size  of  cavities 
formed  at  the  interfaces  increase  as  the  applied  strain  increases  from  1.2%  to  5%,  as  expected. 

Case  2.2:  Inclusion  with  different  radii 
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From  stress  contours  one  observes  several  features'*:  1)  cavities  formed  at  the  interfaces 
between  the  matrix  and  the  four  inclusions  are  not  the  same;  2)  interfacial  failure  process 
at  the  interfaces  are  significantly  different  between  case  2.1  and  case  2.2:  in  case  2.1,  voids 
formed  at  interfaces  will  grow  monotonically,  in  case  2.2  voids  formed  at  interfaces  can  grow 
as  well  as  shrink  due  to  inclusion  interactions.  The  second  observation  is  supported  also  by 
the  difference  in  the  global  stress  strain  responses  for  the  two  subcases,  as  shown  in  Figure 
3.  In  Figure  3,  we  see  clearly  that  the  stress  strain  curves  for  the  two  subcases  are  essentially 
the  same  when  global  strain  is  smaller  than  2.5%.  However,  when  this  strain  exceeds  2.5%, 
kinks  in  these  two  stress  strain  curves  occur  at  different  strain  levels  and  their  profiles  have 
distinct  features.  (Kinks  in  stress  strain  curves  indicate  substantial  debonding  at  interfaces). 

So  when  the  volume  fraction  of  inclusions  is  at  intermediate  level,  interactions  between 
inclusions  are  moderate,  but  effects  of  these  interactions  on  the  damage  evolution  are  easily 
identifiable  from  global  stress  strain  curve. 

High  volume  fraction  of  inclusions  (case  3) 

When  the  volume  fraction  is  as  high  as  68.4%,  one  observed  prounced  effects  of  particle 
interaction  and  of  variation  of  particle  sizes. 

Case  3.1:  Inclusion  with  same  radius 

We  report  the  stress  contours  in  Figures  4a,  5a,  we  see  that  (1)  debonding  occurs  at 
interfaces  near  symmetry  plane  perpendicular  to  the  direction  of  loading,  (2)  separations 
grow  with  the  increase  of  global  strain,  (3)  when  the  global  strain  is  around  10%,  the  interface 
cavities  look  similar  to  those  predicted  by  unit  cell  analyses. 

Case  3.2;  Inclusions  with  differnt  radii 

When  the  sizes  of  the  4  inclusions  are  slightly  different,  it  is  observed  (Figure  4b,  5b  ) 
that  (1)  debonding  occurs  at  the  larger  inclusions  first,  (2)  cavities  formed  at  interfaces  grow 
or  shrink  with  the  increasing  global  strain  (this  is  impossible  in  a  unit  cell  analysis).  This 
second  observation  is  clearly  a  strong  indication  of  the  significant  effect  of  the  interactions 
between  particles  and  also  of  the  significant  effect  of  even  a  slight  variation  in  particle  sizes. 

The  interaction  of  particles  and  variation  of  inclusion  sizes  jointly  control  the  damage 

‘‘To  save  space,  we  do  not  display  stress  contours  here.  More  detailed  numerical  results  are  given  for  case 
3  where  significant  effects  of  size  variation  are  observed 
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evolution  process  as  well  as  the  global  material  response.  The  global  stress  strain  curves 
shown  in  Figure  6  indicates  that  when  significant  debonding  occurs  (indicated  by  the  kink 
at  0.2%  global  strain),  the  material  responses  differ  substantially  between  a  composite  with 
the  same  size  of  inclusions  and  a  composite  possesing  inclusions  of  different  sizes. 

A  16-inclusion  problems  (case  4) 

A  square  with  16  inclusions  was  also  considered  for  which  the  inclusion  sizes  fluctuate 
around  a  mean  of  0.4  rendering  a  volume  fraction  of  50%.  Again,  normal  displacements  are 
described  at  two  opposite  sides  of  the  domain  with  the  remaining  two  sides  being  traction 
free.  For  this  case  the  computation  reaches  4%  global  strain  with  the  largest  normal  strain 
in  the  interior  being  23%;  thereafter  the  algorithm  diverges. 

A  detailed  examination  of  the  results  for  this  case  shows  that  (1)  interaction  between 
particles  control  the  damage  evolution  process;  (2)  a  small  change  in  particles  sizes  has 
a  significant  effect  on  damage  evolution  and  material  response  (see  Figure  6ab  in  Zhong 
and  Knauss  (1997b)).  In  addition,  discrete  segmental  increases  are  observed  in  the  global 
stress-strain  curve. 

The  4  cases  considered  here  indicate  that  the  interaction  of  particles  and  even  the  slight 
variation  in  the  particle  size  control  damage  evolution  in  a  "particulate”  composite  and  its 
local  stress  strain  response. 

4.2  Interaction  of  inclusions:  large  variation  in  inclusion  size 

Again  we  consider  a  square  domain  with  16  inclusions  with  different  radii  as  indicated  in 
Figure  7  (considering  symmetry,  only  a  quarter  of  the  domain  is  displayed).  The  particle 
sizes  are  arranged  as  shown  on  purpose  so  that  one  can  take  advantage  of  the  symmetry 
conditions  involved. 

From  Figures  8  -  9,  we  observe  that  (1)  debonding  occurs  at  interfaces  between  the  largest 
particle  and  matrix  first,  and  subsequently  at  interfaces  of  smaller  particles  (of  descending 
size)  and  matrix;  (2)  there  is  essentially  no  debonding  at  the  interfaces  between  the  smallest 
particles  (about  half  the  size  of  the  large  particle)  and  the  matrix  even  when  the  global  strain 
is  8%  or  so. 

Again,  we  observed  that  the  interaction  between  particles  plays  important  role  in  the 
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damage  (debonding)  evolution  process.  It  seems  that  a  bimodal  distribution  of  particles  in 
an  elastomer  produce  some  distinct  features  of  damage  evolution  and  material  response  (See 
Figure  10,  and  also  compare  Figure  10  with  Figure  7  and  Figure  6a  of  Zhong  and  Knauss 
(1997b)) 

The  strong  effect  of  interaction  between  particles  and  the  significant  effect  of  the  slight 
variation  of  particle  sizes  (for  high  volume  fraction  of  particles)  cannot  be  predicted  by  unit 
cell  analyses  or  damage-incorporating  constitutive  models.  These  effects  might  contribute 
to  failure  mechanisms  in  particulate  composites  with  high  volume  fraction  of  particles  such 
as  solid  propellants  (volume  fraction  70  -  90  %).  Nakamura  et.  al.  (1992)  investigated 
experimentally  the  effect  of  the  particle  size  on  the  fracture  toughness  of  epoxy  resin  filled 
with  spherical  silica.  They  found  that  both  critical  stress  intensity  factor  and  critical  energy 
release  rate  increased  with  (mean)  particle  size.  This  can  be  eeisily  explained  by  the  fact 
that  when  the  (mean)  particle  size  is  large,  significant  damage  (in  the  form  of  interfacial 
debonding)  occurs  at  a  crack  tip.  This  damage  zone  can  consume  stored  strain  energy  at 
the  crack  tip  and  thus  hamper  crack  propagation.  Our  results  seem  to  be  in  line  with  their 
observations  and  support  the  above  explanation  in  the  following  way;  The  proposed  model 
predicted  that  the  larger  the  particle  size  is,  the  easier  is  the  debonding.  The  interfacial 
failure  process  consumes  part  of  the  strain  energy.  This  correlation  between  numerical  and 
experimental  observations  is  indirect  because  we  did  not  consider  the  effect  of  interfacial 
failure  on  fracture  toughness  in  our  numerical  study  yet.  In  addition,  our  results  also  sug¬ 
gested  that  the  standard  deviation  of  particle  size  distribution  should  have  significant  effect 
on  fracture  toughness.  Comparing  global  stress  strain  response  of  composite  of  same  size 
inclusions  and  that  of  composite  with  inclusions  of  slightly  different  size,  the  interfacial  fail¬ 
ure  process  in  the  latter  composite  consumed  more  energy  than  that  in  the  first  composite. 
We  expect  that  the  fracture  toughness  for  the  latter  composite  (small  variation  in  inclusion 
sizes)  will  be  larger. 
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5  Conclusions 


By  characterizing  the  particles,  matrix  and  interfaces  in  a  particulate-filled  elastomer  individ¬ 
ually,  the  whole  failure  process  (interfacial  debonding,  growth  of  voids  formed  at  interfaces) 
is  analyzed  by  the  FEM. 

It  is  shown  that  (1)  When  the  volume  fraction  is  high,  even  small  changes  in  particle  size 
can  have  a  large  effect  on  local  material  response  and  on  damage  evolution  in  particular;  (2) 
When  there  is  a  large  difference  in  particle  sizes(such  as  bimodal  distribution),  separation 
damage  occurs  at  interfaces  between  the  large  particles  and  matrix,  little  or  no  damage 
occurs  at  small  ones.  It  seems  that  bimodal  distribution  of  particles  in  an  elastomer  produce 
some  distinct  features  of  damage  evolution  and  local  material  response. 

Although  the  effects  of  inclusion  size  (small)  variation  on  local  material  response  may 
be  smoothed  out  for  an  ensemble  with  large  number  of  inclusions,  it  seems  likely  that  the 
change  in  local  material  response  will  affect  failure  process  (e.g.  interfacial  debonding)  as 
well  as  continuum  material  properties  such  as  fracture  toughness. 

In  order  to  analyze  the  interaction  of  particles  with  a  macroscopic  crack,  such  as  prop¬ 
agation  of  a  macroscopic  crack  via  the  growth  and  coalescence  of  voids,  the  computational 
approach  must  be  made  more  robust  and  stable.  This  type  of  investigation  may  lead  to  a 
better  understanding  of  failure  mechanism  in  filled  elastomers. 
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Figure  I.  A  discretely  failing  region  embedded 
in  a  linear / nonlinear  material 


Figure  3:  Average  nominal  stress  strain  curves  for  intermediate  inclusion 
volume  fraction. 


Figure  6:  Averge  nominal  stress  strain  curves  for  large  volume  fraction  cases. 
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Contour  of  Sll  for  bimodal  size  distribution  case  (2.2%  global  strain). 
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Figure  10:  Avergare  stress  strain  curve  for  bimodal  size  distribution  case. 


